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This  report  presents  the  results  of  a  three-year  theoretical  and  experimental 
investigation  of  the  response  of  multiphase  porous  media  to  strong  dynamic  loadings. 

A  material  element  model  was  developed  to  automate  the  determination  of  parameter 
values  for  multiphase  material  mdoels  and  to  predict  the  response  of  saturated 
porous  materials  to  arbitrary  stress,  strain  and  fluid  flow  boundary  conditions. 
Laboratory  tests  to  measure  fluid  flow  through  porous  media  were  conducted  over  a 
wide  range  of  pressure  gradients  in  granular  materials  and  porous  rock.  The  results 
of  those  tests  were  used  to  develop  a  unified  fluid  friction  model.  A  laboratory 
apparatus  was  constructed  for  the  purpose  of  studying  wave  propagation  through 
saturated  porous  media.  A  suite  of  dynamic  tests  was  performed  on  a  porous  stainless 
steel  specimen,  and  the  results  of  those  tests  were  compared  with  detailed  multiphase 
numerical  simulations.  The  simulations  were  in  excellent  agreement  with  the 
laboratory  measurements.  With  guidance  from  simulations,  it  was  possible  to  identify 
in  the  test  data,  the  wve  of  the  second  kind.  This  disturbance,  which  propagates 
at  a  speed  significantly  lower  than  the  compresslonai  wave,  appears  to  be  associated 
with  a  surge  of  pore  fluid  within  the  porous  skeleton,  and  has  previously  only  been 
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EXECUTIVE  SUMMARY 


For  more  than  ten  years,  the  Air  Force  Office  of  Scientific  Research  (AFOSR)  has 
supported  research  in  the  field  of  multiphase  mechanics,  i.e.  the  mechanical  response  of 
materials  consisting  of  a  porous  skeleton  of  solid-phase  material  having  its  pore  space  filled  with 
liquid  or  a  mixture  of  liquid  and  gas.  For  the  Department  of  Defense  (,DoD),  the  primary 
motivation  for  such  efforts  is  the  study  of  the  effects  of  explosions  in  fully  or  partially  saturated 
soils  and  porous  rocks,  which  exhibit  significantly  different  behavior  than  equivalent  dry 
materials.  This  fundamental  research  in  multiphase  mechanics  is  also  applicable  to  other 
problems  such  as  earthquake  induced  liquefaction  and  shock-wave  therapies  for  kidney  stones. 
The  results  of  ^he  past  decade  of  AFOSR-sponsored  research  has  contributed  significantly  to 
numerous  DoD  studies,  including: 

•  an  evaluation  of  nuclear  craters  in  the  Pacific  which  led  to  a  revision  in  the 
perceived  cratering  effectiveness  of  nuclear  weapons; 

•  evaluations  of  the  vulnerability  of  deep  underground  structures  in  saturated  porous 
rock  geologies; 

•  development  of  simplified  weapons  effects  prediction  techniques  for  use  by 
targeting  analysts  throughout  the  DoD  community; 

•  analysis  of  the  vulnerability  of  U.S.  Air  Force  aircraft  shelters  and  other 
hardened  facilities  sited  in  saturated  soils  to  attack  from  conventional  weapons; 

•  development  of  countermeasures  to  enhance  the  survivability  of  hardened 
structures  in/on  saturated  soil;  and 

•  analysis  and  experimental  support  for  the  U.S.  Navy’s  surf  zone  mine 
countermeasures  research. 

Analysis  of  the  response  of  saturated  porous  materials  is  complicated  by  the  fact  that  they 
are  composed  of  multiple  constituents,  e.g.  sand  and  water,  having  fundamentally  different 
mechanical  response  characteristics.  A  computation  of  the  deformation  of  the  multiphase 
material  under  load  must  account  for  the  partitioning  of  load  between  the  solid  and  fluid  phases. 
In  addition  to  the  deformation  of  the  solid  skeleton  caused  by  the  load  it  carries,  it  is  also 
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deformed  by  the  pore  fluid  pressure  acting  on  it.  Thus,  a  theoretically  rigorous  formulation 
must  account  for  the  deformations  of  both  materials  and  the  coupling  between  them  while 
maintaining  full  compatibility  of  deformations  as  well  as  tracking  mass  transfer  of  the  fluid 
phase.  All  geologic  and  cemenlitions  materials  exhibit  frictional  strength  characteristics,  i.e. 
over  a  wide  stress  range,  their  shear  strength  increases  with  increasing  normal  stress.  However, 
in  the  case  of  fluid-saturated  materials,  only  that  component  of  normal  stress  that  is  carried  by 
the  solid  skeleton  is  effective  in  increasing  its  shear  strength.  This  "effective  stress"  concept, 
which  was  developed  by  Terzaghi  over  50  years  ago,  was  largely  ignored  by  the  DoD  weapons 
effects  community  until  the  results  of  the  AFOSR  research  began  to  influence  the  community 
in  the  last  few  years. 


The  research  effort  described  in  this  report,  which  was  conducted  by  Applied  Research 
Associates,  Inc.  (ARA),  was  designed  to  further  advance  the  state  of  knowledge  in  multiphase 
mechanics.  It  was  a  rather  ambitions  program  with  three  major  objectives: 


1.  development  of  a  Material  Element  Model  (MEM)  to  automate  the  determination 
of  parameter  values  for  multiphase  material  models  and  to  predict  the  response 
of  saturated  undrained  porous  materials  to  specified  stress  and  strain  inputs; 

2.  verification  and  calibration  of  a  fluid  friction  model  for  soils  and  rocks  covering 
a  wide  range  of  pore  pressure  gradients,  and  correlation  of  the  fluid  friction 
parameters  with  the  microscopic  structure  of  the  pore  space;  and 

3.  development  of  a  better  understanding  of  the  fundamental  mechanics  of  wave 
propagation  in  saturated  porous  media  through  theoretical,  numerical  and 
experimental  studies  of  wave  propagation  in  confined  saturated  porous  bars. 


While  the  principles  of  multiphase  mechanics  and  the  effective  stress  law  have  been 
known  for  many  years,  their  application  to  problems  of  interest  to  the  DoD  community  and  the 
accomplishment  of  further  research  were  impeded  by  the  lack  of  analysis  tools.  Even  the 
simplest  problems  in  multiphase  mechanics  are  analytically  tractable  only  with  numerical 
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methods  on  computers.  Because  of  the  direct  influence  of  pore  fluid  pressure  on  strength  and 
the  sensitivity  of  pore  pressure  to  changes  in  volume,  it  is  extremely  important  in  multiphase 
calculations  that  the  volume  change  behavior  of  the  solid  skeleton  under  shear  and  normal 
loading  be  modeled  correctly.  Under  this  contract,  ARA  developed  a  Material  Element  Model 
(MEM)  to  support  constitutive  modeling  of  multiphase  materials.  It  is  a  quasi-static  finite 
element  program  that  was  specifically  designed  to  facilitate  implementation  of  new  skeleton 
constitutive  models.  It  provides  a  very  convenient  platform  for  constitutive  model  parameter 
evaluation  and  parameter  sensitivity  studies.  The  resulting  material  model  parameter  values  are 
directly  transportable  to  dynamic  analysis  programs,  a  highly  desirable  feature  in  view  of  the 
fact  that  the  modern  elastoplastic  constitutive  models  may  have  as  many  as  40  parameters. 
MEM  includes  an  extremely  robust  set  of  boundary  conditions,  making  it  possible  to  exercise 
material  models  over  virtually  any  combination  stress,  strain,  and/or  fluid  flow  conditions. 

A  description  of  the  formulations  implemented  in  MEM  is  presented  in  Section  2.  An 
advanced  three-invariant  elastoplastic  constitutive  model  which  is  included  in  MEM  for  modeling 
of  multiphase  materials  is  described  in  Section  3,  along  with  a  discussion  of  the  procedures  for 
fitting  the  model  parameters  for  a  granular  material.  An  important  application  of  MEM, 
presented  in  Section  5,  involves  the  analysis  of  a  conventional  triaxial  compression  test  on  a 
porous  limestone  specimen  in  the  saturated  undrained  condition.  When  the  test  data  were 
interpreted  using  conventional  approaches,  the  results  led  to  the  apparent  paradox  that  the  pore 
fluid  pressure  in  the  specimen  remained  constant  while  the  volume  of  the  specimen  increased 
under  triaxial  compression  loading  at  constant  confining  pressure.  A  study  involving  carefully 
conducted  tests  and  a  detailed  analysis  using  MEM  led  to  the  conclusion  that  the  conventional 
triaxial  compression  test  applied  to  saturated  undrained  material  is  fatally  flawed,  and  that  even 
for  dry  materials,  it  can  lead  to  erroneous  measures  of  volume  change  under  shear  loading.  It 
has  long  been  recognized  that  triaxial  compression  test  specimens  deform  nonuniformly  owing 
to  the  additional  effective  radial  confinement  resulting  from  end  friction.  However,  the 
traditional  approach  has  been  to  assume  that  error  occurs  only  at  the  ends  of  the  specimen,  and 
that  the  behavior  of  the  center  third  of  the  specimen  is  representative  of  the  material’s  response 
to  the  idealized  triaxial  compression  loading.  In  fact,  both  the  test  and  the  MEM  analysis 
showed  that  the  deformation  of  the  specimen  was  nonuniform  in  both  the  axial  and  radial 
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directions,  with  the  large  expansive  radial  deformation  in  the  center  of  the  specimen 
accompanied  by  a  larger  than  average  compressive  axial  strain.  Thus,  the  traditional  approach 
of  volume  strain  computation  from  the  radial  deformation  at  the  center  of  the  specimen  and  the 
axial  deformation  averaged  over  the  entire  length  substantially  overstates  the  volumetric 
expansion  of  the  specimen.  In  subsequent  work,  motivated  by  the  results  of  this  study,  a 
procedure  was  developed  to  reduce  end  friction  using  lubricating  membranes.  Using  that  test 
technique,  we  have  demonstrated  that,  at  confining  pressures  above  the  brittle-ductile  transition, 
porous  limestone  will  undergo  axial  strains  up  to  15%  with  no  tendency  to  dilate.  This  contrasts 
with  the  conventional  approach  which  inferred  dilation  as  a  result  of  the  nonuniform  strain  fields 
in  the  test  specimen.  As  a  result  of  this  study  of  test  techniques  for  saturated  porous  materials 
using  MEM,  modifications  have  been  implemented  to  rock  testing  techniques  in  laboratories  at 
both  ARA  and  WES.  These  results  were  also  presented  at  an  international  conference  on 
constitutive  modeling  (Chitty,  et  al.,  1991). 

The  second  major  thrust  of  this  research  was  in  the  area  of  fluid  flow  through  granular 
materials  and  porous  rocks  at  high  pressure  gradients.  Dynamic  loading  of  a  saturated  porous 
material  will  induce  pore  pressure  gradients.  The  development  of  pore  pressure,  and  hence  the 
strength  of  the  saturated  material,  under  these  conditions  is  strongly  influenced  by  its  fluid  flow 
characteristics.  Therefor,  to  properly  model  such  dynamic  loading,  it  is  necessary  to  analytically 
treat  the  fluid  flow  as  well  as  the  compressibility  and  strength  of  the  material.  We  have 
proposed  and  validated  an  equation  of  motion  for  fluid  flow  that  includes  inertial  terms  as  well 
as  flow  resistance  terms  that  are  linear  and  quadratic  in  flow  velocity.  While  the  drag  equation 
involving  the  linear  and  quadratic  terms  was  first  proposed  by  Forchheimer  near  the  beginning 
of  this  century  and  has  been  used  in  chemical  engineering,  it  was  virtually  unknown  in  the  civil 
engineering  literature  and  had  not  been  validated  at  high  pressure  gradients.  We  initially 
believed  that  the  appearance  of  the  quadratic  term  was  related  to  turbulent  flow,  analogous  to 
flow  in  pipes.  However,  analysis  of  test  data  for  both  granular  materials  and  porous  rocks  has 
shown  that,  unlike  pipes  in  which  there  is  a  change  from  linear  to  quadratic  dependence  on 
velocity  at  the  transition  from  laminar  to  turbulent  flow,  both  terms  are  present  over  the  entire 
range  of  pressure  gradients  in  porous  materials.  We  now  believe  that  the  linear  term  is  a  result 
of  drag  forces  related  to  the  size  of  the  narrow  passages  or  throats  through  which  the  fluid  must 


VI 


flow,  and  the  quadratic  term  is  a  result  of  changes  in  momentum  due  to  direction  changes  as  the 
fluid  flows  along  a  tortuous  path  through  the  material.  The  Darcy  equation  in  conventional  soil 
mechanics  is  equivalent  to  retaining  the  linear  term  and  neglecting  the  quadratic  term.  This 
approximation  has  served  well  over  the  years  because  the  effect  of  the  quadratic  term  becomes 
small  in  comparison  with  the  linear  at  the  low  velocities  (pressure  gradients)  encountered  in 
hydrology  and  in  conventional  civil  engineering  applications.  We  have  adopted  a  form  of  the 
flow  equation  in  which  the  flow  coefficients  have  dimensions  involving  only  powers  of  length 
and  are  functions  only  of  the  properties  of  the  porous  material,  independent  of  fluid  properties. 
The  fluid  density  and  viscosity  appear  explicitly  in  the  equation.  Flow  tests  were  performed  on 
granular  materials  with  a  range  of  sizes  and  shapes  and,  under  sponsorship  of  the  Defense 
Nuclear  Agency  (DNA),  on  porous  limestone.  Relationships  have  been  developed  to  relate  the 
fluid  flow  coefficients  to  the  geometry  of  the  porous  materials.  The  linear  flow  coefficient  is 
related  to  representative  throat  dimensions  of  the  porous  material  by  a  dimensionless  constant 
over  the  whole  range  of  granular  materials  and  rock  tested.  A  similar  relationship  was 
developed  for  the  quadratic  flow  coefficient,  but  it  is  found  to  be  applicable  only  to  the  granular 
materials.  Since  the  tortuosity  in  rock  is  fundamentally  different  than  in  granular  materials,  this 
is  consistent  with  the  view  that  the  quadratic  coefficient  is  related  to  the  inertial  forces  associated 
with  direction  changes.  The  fluid  flow  test  results  and  theoretical  development  of  the  flow 
model  are  documented  in  Section  4  of  the  report. 

The  final  major  objective  of  this  research  was  to  perform  laboratory  wave  propagation 
tests  on  saturated  porous  materials  for  the  purpose  of  validating  the  analytical  and  numerical 
formulations,  and  observing  their  fundamental  mechanical  response  characteristics.  The  early 
part  of  the  research  effort  involved  parallel  efforts  to  design  the  laboratory  apparatus  and 
perform  numerical  analyses  of  the  planned  experiments.  The  calculational  portion  of  the  effort 
was  performed  using  the  MultiPhase  Dynamic  Analysis  Program  (MPDAP)  which  was  developed 
under  a  previous  AFOSR-sponsored  effort.  For  this  work,  MPDAP  which  was  originally  an 
implicit  finite  element  program  was  converted  to  have  an  explicit  integration  option  which  will 
more  efficiently  treat  wave  propagation  problems,  as  described  in  Section  2.  Our  approach  to 
design  of  the  test  apparatus  was  similar  in  appearance  to  a  Split  Hopkinson  Pressure  Bar 
(SHPB).  A  striker  bar  is  fired  from  a  gas  gun  and  impacts  the  incident  bar  inducing  a 
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compressive  wave  which  is  propagated  into  the  test  specimen.  In  contrast  to  the  SHPB,  the 
Wave  Propagation  Bar  (WPB)  makes  measurements  directly  on  a  test  specimen  which  is  many 
bar  diameters  in  length.  The  ARA  WPB  uses  a  specimen  2  inches  in  diameter  and  up  to  24 
inches  long.  Because  of  the  necessity  to  saturate  the  test  specimen,  it  was  originally  proposed 
to  confine  the  specimen  within  pressurized  gas  in  a  concentric  pressure  vessel,  much  like  a 
triaxial  compression  apparatus.  However,  the  MPDAP  numerical  calculations  revealed  that  the 
experiment  as  planned  had  significant  shortcomings.  Because  the  confining  pressure  remained 
essentially  constant  throughout  the  test,  the  pore  pressure  generated  by  the  propagating  wave  was 
prematurely  relieved  by  a  relief  wave  propagating  in  from  the  surface  of  the  porous  specimen. 
This  discovery  motivated  further  calculations  to  find  a  solution  to  the  problem  and  caused  a 
delay  in  construction  of  the  test  equipment.  Eventually,  based  on  the  calculatior.al  results,  an 
apparatus  was  designed  and  constructed  in  which  the  confining  vessel  consisted  of  a  thick 
(approx.  2  inches)  steel  sleeve  with  a  thin  (approx.  0.2  inch)  liquid  filled  gap  between  its  inner 
wall  and  the  test  specimen.  A  porous  metal  bar  was  chosen  as  a  test  specimen  because  it 
appeared  to  have  a  combination  of  skeleton  response  and  fluid  flow  properties  which  would 
accentuate  the  influence  of  relative  fluid  motion  on  the  wave  propagation  characteristics.  The 
design  of  the  equipment  and  the  preliminary  calculations  are  presented  in  Section  6. 

The  WPB  test  apparatus  was  constructed  and  a  series  of  tests  was  run  on  a  specimen  of 
porous  stainless  steel  in  both  the  saturated  and  dry  states,  as  described  in  Section  7. 
Measurements  were  made  of  dynamic  skeleton  strains  at  three  locations  on  the  specimen  and  of 
dynamic  pore  pressure  respors?  at  one  location.  The  experimental  apparatus  produced  highly 
repeatable  measurements  of  excellent  quality.  As  also  documented  in  Section  7,  the  test  results 
were  analyzed  using  MPDAP  calculations  which  simulated,  as  closely  as  possible,  the  lab  tests. 
There  is  excellent  qualitative  and  quantitative  agreement  between  the  calculations  and  the 
laboratory  measurements.  The  calculational  results  clearly  exhibit  a  secondary  disturbance 
propagating  at  a  relatively  slow  speed  behind  the  compressional  wave.  This  wave  of  the  second 
kind  is  predicted  by  analysis  to  exist  in  saturated  porous  materials  but  has  not,  to  our  knowledge, 
previously  been  observed.  Guided  by  the  calculations,  we  were  able  to  identify  this 
phenomenon,  which  appears  to  be  related  to  a  surge  of  pore  fluid  within  the  porous  skeleton, 
in  the  laboratory  test  data.  Thus,  we  believe  that  we  have  observed  a  phenomenon  in  the 
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laboratory  which  has  previously  only  been  predicted  by  analysis.  The  excellent  overall 
agreement  between  experimental  wave  propagation  test  data  and  the  numerical  calculations  of 
these  experiments  provides  strong  validation  of  these  analytical  techniques.  It  should  provide 
confidence  in  their  application  to  the  much  more  complex  real-world  problems  of  interest  to  the 
DoD  and  others. 
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SECTION  1 
INTRODUCTION 


This  report  presents  the  results  of  a  research  effort  directed  at  improving  the  fundamental 
understanding  of  the  dynamic  response  of  fully  or  partially  fluid-saturated  porous  materials.  The 
work  was  performed  by  Applied  Research  Associates,  Inc.  (ARA)  under  Contract  No.  F49620- 
89-C-0037  with  the  Air  Force  Office  of  Scientific  Research  (AFOSR).  The  laboratory  research 
was  conducted  by  the  Materials  Testing  Laboratory  of  ARA’s  New  England  Division  in  South 
Royalton,  Vermont,  which  was  also  responsible  for  the  data  analysis  and  numerical  simulations. 
An  ARA  subcontractor,  COMTEC  Research  of  Clifton,  Virginia  developed  the  multiphase 
analysis  programs,  MEM  and  MPDAP. 

1.1  BACKGROUND 

For  more  than  ten  years,  AFOSR  has  supported  a  number  of  researchers  in  the  field  of 
multiphase  mechanics,  with  particular  emphasis  on  blast  induced  liquetaction.  This  research  has 
significantly  advanced  the  state  of  knowledge  in  this  field  of  geomechanics  and  has  provided 
strong  support  for  Air  Force  and  Department  of  Defense  studies  of  weapons  effects  in  saturated 
geologies.  Studies  which  this  work  has  supported  include: 

•  an  evaluation  of  nuclear  craters  in  the  Pacific  which  led  to  a  revision  in  the 
perceived  cratering  effectiveness  of  nuclear  weapons; 

•  evaluations  of  the  vulnerability  of  deep  underground  structures  in  saturated  porous 
rock  geologies; 

•  development  of  simplified  weapons  effects  prediction  techniques  for  use  by 
targeting  analysts  throughout  the  DoD  community; 
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•  analysis  of  the  vulnerability  of  U.S.  Air  Force  aircraft  shelters  and  other 
hardened  facilities  sited  in  saturated  soils  to  attack  from  conventional  weapons; 

•  development  of  countermeasures  to  enhance  the  survivability  of  hardened 
structures  in/on  saturated  soil;  and 

•  analysis  and  experimental  support  for  the  U.S.  Navy’s  surf  zone  mine 
countermeasures  research. 


1.2  OBJECTIVES 

The  research  described  in  this  report  was  designed  to  further  advance  the  state  of 
knowledge  in  the  are  of  dynamic  multiphase  mechanics.  There  are  three  major  objectives  of  this 
research: 

1.  development  of  a  Material  Element  Model  (MEM)  to  automate  the  determination 
of  parameter  values  for  multiphase  material  models  and  to  predict  the  response 
of  saturated  undrained  porous  materials  to  specified  stress  and  strain  inputs; 

2.  verification  and  calibration  of  a  fluid  friction  model  for  soils  and  rocks  covering 
a  wide  range  of  pore  pressure  gradients,  and  correlation  of  the  fluid  friction 
parameters  with  the  microscopic  structure  of  the  pore  space;  and 

3.  development  of  a  better  understanding  of  wave  propagation  mechanics  in 
saturated  porous  media  through  theoretical,  numerical  and  experimental  studies 
of  wave  propagation  in  confined  saturated  porous  bars. 


1.3  SCOPE 


Model  development  work  under  this  effort  is  summarized  in  Sections  2  and  3.  The 
formulation  of  the  computer  program,  MEM,  a  tool  for  the  development  of  multiphase 
constitutive  models  for  fully  and  partially  saturated  porous  materials,  and  the  extension  of 
MPDAP,  a  finite  element  program  for  multiphase  dynamic  analysis,  to  allow  for  explicit 
integration,  making  it  more  efficient  for  wave  propagation  problems  are  detailed.  The 
formulations  of  the  multiphase  mechanics  are  presented  in  Section  2.  Section  3  describes  the 
constitutive  relationships  used  to  model  the  individual  constituents  in  a  multiphase  problem, 
including  the  porous  skeleton,  the  solid  grains,  and  the  pore  fluid.  The  source  codes  and  user 
guides  for  MEM  and  MPDAP  were  presented  in  the  first  two  Annual  Technical  Reports 
submitted  under  this  contract  (Blouin,  et.  al;  1990,  1991) 

Laboratory  tests  to  support  development  of  fluid  flow  models  for  porous  media  were 
conducted  at  a  wide  range  of  pressure  gradients  in  granular  materials  with  a  range  of  grain  sizes 
and  in  porous  rock.  The  results  of  fluid  flow  research  are  presented  in  Section  4  along  with 
analysis  which  results  in  a  unified  fluid  friction  model.  Section  5  presents  the  results  of  an 
analysis  of  the  undrained  triaxial  compression  test  on  saturated  undrained  triaxial  compression 
test  on  a  porous  limestone  specimen  using  data  and  models  from  the  previous  sections.  The  test 
results  are  compared  with  a  detailed  multiphase  numerical  simulation  of  the  test.  In  addition  to 
exercising  and  validating  the  numerical  modeling  techniques,  this  analysis  identified  fundamental 
problems  with  the  test  and  provided  suggestions  for  correcting  them. 

In  order  to  further  study  the  dynamic  respon.se  of  saturated  porous  materials  and  validate 
procedures  for  numerical  simulation  of  that  response,  a  laboratory  apparatus  was  designed  and 
fabricated  to  study  wave  propagation  through  fluid-confined  saturated  porous  materials.  Section 
6  documents  the  development  of  that  apparatus,  including  a  suite  of  supporting  numerical 
simulations.  A  series  of  tests  was  performed  on  a  bar  of  porous  (sintered)  stainless  steel  in  both 
the  dry  and  fully  saturated  conditions.  Section  7  describes  those  tests  and  presents  the  test 
results  along  with  comparisons  to  a  numerical  simulation  of  one  of  the  tests  that  was  performed 
using  the  models  described  in  the  previous  sections. 


SECTION  2 

NUMERICAL  MODELING  OF  MULTI-PHASE  MEDIA 


The  differential  equadons  governing  wave  propagation  in  saturated  porous  media  are  so 
complex  that  closed-form  analytical  solutions  are  totally  impractical  for  all  but  the  simplest  one¬ 
dimensional  geometries.  The  analytical  difficulty  is  further  exacerbated  in  the  case  of  nonlinear 
fluid  or  skeleton  properties,  a  great  variety  of  numerical  algorithms  have  been  developed  to 
model  the  mechanical  deformation  of  materials  under  static  and  dynamic  loading.  However, 
numerical  techniques  that  can  treat  saturated  porous  media  have  not  been  well  developed.  This 
section  describes  two  different  numerical  procedures  for  analysis  of  saturated  or  partially 
saturated  porous  materials. 

The  Material  Element  Model  (MEM)  is  a  quasi-static  finite  element  program  that  was 
developed  under  this  research  effort  to  facilitate  constitutive  model  development  for  saturated 
porous  materials.  It  models  deformation  and  fluid  flow  under  arbitrary  loading,  specified  in 
terms  of  either  stress  or  strain,  and  it  is  designed  to  simplify  inclusion  of  new  skeleton 
constitutive  models.  Because  pore  pressure,  and  hence  strength,  are  strongly  dependent  on  the 
volumetric  behavior  of  the  porous  skeleton,  accurate  constitutive  modeling  of  the  skeleton  is 
particularly  important  in  the  analysis  of  saturated  materials. 

Multi-phase  analyses  require  greater  computation  time  and  computer  memory  than  single¬ 
phase  dynamic  analyses  because  of  the  required  additional  degrees  of  freedom  to  represent  the 
motion  of  pore  fluid.  Consequently,  computational  efficiencies  are  an  essential  consideration 
in  the  formulation  of  solution  schemes  for  real  problems  of  interest  in  weapons  effects  and 
earthquake  engineering.  This  section  describes  an  efficient  explicit  formulation  for  the  multi¬ 
phase  dynamic  analysis.  In  contrast  to  implicit  solution  techniques,  the  explicit  formulation  does 
not  require  assembly  and  solution  of  a  global  stiffness  matrix.  For  numerical  stability,  the 
explicit  approach  requires  time  steps  that  are  no  greater  than  the  wave  transit  time  of  the 
smallest  elements.  While  this  is  smaller  than  the  step  size  required  for  numerical  stability  in  an 
implicit  calculation,  wave  propagation  problems  typically  require  time  steps  of  that  size  to 
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capture  the  physics  of  interest.  Thus,  the  explicit  solution,  with  its  greatly  reduced  numerical 
complexity,  is  the  preferred  method  for  real  wave  propagation  problems. 

2.1  MATERIAL  ELEMENT  MODEL 

MEM  is  a  finite  element  program  that  models  the  quasi-static  response  of  saturated 
porous  materials  to  prescribed  stress,  strain  and  pore  pressure  boundary  conditions.  As  currently 
implemented,  MEM  has  both  single  element  and  multiple  element  options.  The  single  element 
option  includes  full  three  dimensional  capability.  The  multi  element  option  is  two-  dimensional, 
and  both  options  model  pore  fluid  flow  between  elements.  MEM  is  written  in  updated 
Lagrangian  format  to  enable  modeling  of  nonlinear  geometric  distortions  including  changes  in 
porosity. 

MEM  is  used  to  compute  the  stress,  strain  and  pore  pressure  response  of  saturated, 
nonlinear  porous  materials.  As  depicted  in  Figure  2.1,  MEM  formulates  the  constitutive 
relaitonships  for  9  saturated  porous  materials  in  terms  of  the  nonlinear  response  of  the  individual 
material  components,  including  the  stress-strain  response  of  the  porous  skeleton  in  the  drained 
or  unsaturated  condition,  and  the  pressure-volume  relationships  for  the  solid  grains  and  the  pore 
fluid.  These  are  combined  using  fully  coupled  two-phase  theory  as  described  by  Kim,  Blouin, 
Chitty  and  Merkle  (1988)  and  Simons  (1989).  In  the  fully  coupled  model,  the  overall  material 
response  is  governed  by  the  response  of  the  material  skeleton  to  the  applied  effective  stresses 
as  well  as  the  response  of  the  skeleton  to  compression  of  the  solid  grains  by  the  effective  stresses 
and  pore  fluid  pressure.  In  addition  to  the  material  property  inputs,  a  set  of  boundary  conditions 
and/or  loading  histories  is  input. 

MEM  is  highly  versatile  in  that  it  can  handle  a  wide  variety  of  input  boundary  conditions. 
As  shown  in  Figure  2.2,  individual  stress,  strain  and  pore  pressure  inputs  can  be  prescribed  on 
any  face  of  the  material  element.  In  addition,  boundary  conditions  can  include  prescribed  stress 
followed  by  prescribed  strain,  or  prescribed  strain  followed  by  prescribed  stress.  Pore  fluid 
boundary  conditions  include  drained,  undrained,  and  prescribed  pressure,  as  well  as  drained 
followed  by  undrained  and  vice  versa.  A  good  example  of  a  test  using  a  combination  of  these 


2-2 


boundary  conditions  is  the  shock  consolidation  test  described  by  Kim,  Blouin  and  Timian  (1987). 
In  this  uniaxial  strain  test,  a  saturated  material  undergoes  an  undrained  load-unload  cycle 
followed  by  a  period  of  drained  consolidation  under  constant  effective  stress  and  pore  pressure. 

Two  nonlinear  skeleton  material  models  are  currently  implemented  in  MEM,  a 
generalized  elastic-perfectly  plastic  model  described  by  Kim,  Piepenburg,  and  Merkle  (1985) 
(and  updated  by  Kim  and  Blouin,  1990),  and  the  ARA  Three  Invariant  Conic  Model  described 
by  Merkle  and  Dass  (1985).  The  elastic-plastic  model  requires  seven  material  parameters  to 
completely  describe  the  material  properties.  The  failure  and  yield  surfaces  are  coincident  and 
are  specified  as  one  of  three  options:  Von  Mises,  linear  Mohr-Coulomb,  or  nonlinear  Mohr- 
Coulomb  with  Von  Mises  limit.  The  slope  of  the  nonlinear  failure  envelope  for  rock  is 
prescribed  according  to  Hoek  and  Brown’s  empirical  data  table  (1982)  as  a  function  of  rock  type 
and  in  situ  rock  mass  quality  (Kim,  Blouin  and  Timian,  1986).  Beneath  the  failure  surface  the 
shear  behavior  is  linear-elastic  and  the  pressure-volume  response  is  nonlinear-inelastic.  Along 
the  failure  surface  dilatancy  is  prescribed  according  to  the  associated  flow  rule. 

The  ARA  Three  Invariant  Conic  Model  (ARA  3-1)  is  a  sophisticated  elastic-plastic  model 
with  compressive  and  expansive  work  hardening.  It  requires  about  40  parameters  to  fully 
describe  the  material  model  properties.  These  are  determined  by  computational  fits  to  laboratory 
test  data  as  described  in  Section  3.  Dilatancy  is  controlled  by  a  plastic  potential  function  and 
is  generally  not  associated  with  the  yield  surface.  Details  of  the  current  ARA  3-1  model  are 
described  in  Section  3,  including  improvements  that  have  been  implemented  to  better  describe 
material  behavior. 

The  multi-element  option  of  MEM  permits  the  study  of  two-dimensional  problems  with 
pore  fluid  flow,  e.g.,  two  dimensional  consolidation,  the  simulation  of  saturated  undrained 
laboratory  tests  to  study  the  effects  of  pore  fluid  migration  on  sample  response,  and  simulation 
of  permeability  tests.  The  multi-element  option  has  all  the  boundary  condition  options  listed  for 
the  single  element  in  Figure  2.2  with  the  exception  of  the  changes  from  one  boundary  condition 
to  another. 
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Other  applications  of  MEM  include  the  following: 

•  provide  a  general-purpose  tool  for  testing  engineering  material  models, 

•  provide  a  means  of  constructing  sophisticated  two-phase  material  models  from 
drained  laboratory  test  data  and  of  predicting  equivalent  single-phase  response  to 
prescribed  loadings, 

•  define  the  shear  strength  of  saturated  porous  materials  for  prescribed  undrained 
loadings, 

•  verify  the  two-phase  models  through  comparison  to  laboratory  test  data  obtained 
from  less  conventional  strain  and  stress  paths, 

•  permit  analysis  and  prediction  of  liquefaction  phenomena  under  complex  loading 
conditions, 

•  perform  analysis  of  standard  and  non-standard  laboratory  tests  to  gain  a  better 
understanding  of  test  results  and  to  develop  modifications  to  methods  used  in  the 
determination  of  material  model  parameters. 

2.1.1  Theory 

The  theoretical  formulations  used  in  MEM  are  similar  to  those  used  in  the  multiphase 
dynamic  code  MPDAP  described  by  Kim,  Blouin,  Chitty  and  Merkle  (1988).  The  formulation 
is  briefly  summarized  in  the  following  subsections. 

Full  documents  along  with  MEM  Users  Manual  and  the  listing  of  the  program  are 
presented  by  Blouin,  et  al.,  1990. 

Notation 

Note  that  positive  signs  are  used  for  elongation  and  tension.  A  comma  denotes 
differentiation  with  respect  to  the  subsequent  indices  and  the  superposed  dot  denotes 
diffemtiation  with  respect  to  time.  Prime  indicates  effective  stress  or  effective  pressure. 
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skeleton  displacement 
absolute  fluid  displacement 

apparent  fluid  displacement  relative  to  the  solid  skeleton 

total  stress 

effective  stress 

total  pressure 

effective  pressure 

pore  fluid  pressure 

pore  fluid  pressure  gradient  vector 

skeleton  strain 

element  nodal  skeleton  displacement  vector 

element  nodal  pore  fluid  pressure 

global  nodal  skeleton  displacement  vector 

global  nodal  fluid  pressure 

applied  boundary  traction 

specified  boundary  flow  velocity  (flux) 

body  force  vector  (generally  equals  gravity  force) 

elasto-plastic  stress-strain  matrix  for  skeleton 

unit  vector  (1)^  =  <1  1  1  0  0  0> 

porosity 

pore  fluid  compressibility  =  l/Kf 
compressibility  of  solid  grains 

compressibility  of  soil-water  mixture  with  zero  effective  stress  =  1/K„ 

bulk  modulus  of  soil-water  mixture  with  zero  effective  stress 

bulk  modulus  of  skeleton 

bulk  modulus  of  pore  fluid 

constrained  modulus 

bulk  mass  density  of  mixture 

dry  mass  density  of  skeleton 

fluid  mass  density 

mass  density  of  solid  grains 
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H  :  pore  fluid  viscosity 

a  :  linear  pore  fluid  flow  coefficient  (defined  in  Section  4) 

0  :  quadratic  pore  fluid  flow  coefficient 

[MJ  :  mass  matrix 

[Kj]  ;  tangent  skeleton  stiffness  matrix 

[C]  :  coupling  matrix  between  solid  skeleton  and  pore  fluid 

[E]  ;  pore  fluid  compressibility  matrix 

[H]  :  fluid  friction  energy  dissipation  matrix 

{F}  :  nodal  force  vector 

{R}  :  internal  resistance  force  vector 

{Q}  :  equivalent  boundary  flow  vector 

Formulations 


The  current  version  of  MEM  can  only  handle  the  quasi-static  response  of  the  saturated 
porous  medium.  The  formulations  in  MEM  are  obtained  by  dropping  the  inertia  terms  out  of 
the  general  equations  that  were  derived  for  MPDAP  by  Kim,  Blouin,  Chitty,  and  Merkle  (1988). 


The  first  global  equation  is  based  on  the  equilibrium  of  the  total  stresses  with  the  applied 
boundary  tractions  and  is  given  by 


[ATj  (iiJ,  *  [q  •  {p,l. 


(2.1) 


where: 


[Bf  dv 
[MJ  =  E4  iN\H9-n9j)  [N]  dv 
[JCj]  =  r/,  [Bf  [D^n  IB]  dv 
{f}„  =  E/,  [N]^  {7}„  ds 

[C]  =  E4  [Bf  fill  -  {l}j  <G>  dv 
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The  second  global  equation  is  based  on  the  fact  that  the  internal  fluid  movements  relative 
to  the  solid  skeleton  are  compatible  with  the  specified  boundary  flux  and  is  given  by 


tq’-  (Aa,  »  (-[£]  -  ^  [H]j  (AS.  . 


where; 


[£]  =  EI<G>^  ^  (iF  [D^n  <1} 


g 

[H]  =  E;,  [Af  k^[A]  dv 

iPX  =  -H-i  -  ~  "  ^QK) 

{<?}  =  E  <G>^  Q  ds 

;  /  f  Af  Pf  S  . 
it  =  —  +  |w  . 

\  a  e 
w.  =  niUi  -  li.) 

(1  ■ 


<G>  dv 


-1 


(2.2) 


The  shape  functions  ,  [N]  and  <G> ,  and  their  derivatives  [B]  and  [A],  respectively,  in 
Equations  2.1  and  2.2  relate  the  field  variables  within  each  element  to  .he  element  nodal  values 
as  follows: 

{«}  =  [iV]  (2.3) 

{e}  =  [B] 

It  =  <G>  (ti}^ 

=  [A]  (it}^ 
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2.3  EXPLICIT  FINITE  ELEMENT  FORMULATIONS  FOR  MULTI-PHASE 

DYNAMIC  ANALYSIS 

Kim,  et  al.  (1988)  and  Blouin,  et  al.  (1991)  describe  the  implementation  of  the  multi¬ 
phase  dynamic  equilibrium  equations  using  implicit  integration  schemes.  In  the  implicit  method, 
the  element  stiffness,  damping,  and  mass  matrices  are  assembled  to  form  a  generalized  structural 
stiffness  matrix.  The  linearized  structural  equilibrium  equations  are  then  solved  simultaneously 
at  each  step  in  the  calculation.  For  large  scale  nonlinear  problems,  the  implicit  method  requires 
considerably  longer  computation  at  each  time  step,  mostly  for  solving  linear  equations,  and 
considerably  more  computer  memory  to  store  the  generalized  structural  stiffness  matrix,  than 
a  comparable  problem  solved  using  an  explicit  formulation.  In  many  structural  problems,  the 
penalty  in  computation  time  required  by  the  implicit  approach  at  each  step  is  more  than  offset 
by  a  greatly  reduced  number  of  time  steps  required  to  solve  the  problem.  However,  in  the  wave 
propagation  problems  of  interest  to  members  of  the  DOD  community,  time  steps  on  the  order 
of  the  wave  transit  time  across  an  element  are  required  to  simulate  the  physical  processes  with 
adequate  fidelity.  In  this,  the  explicit  approach  which  requires  approximately  the  same  number 
of  time  steps,  has  an  advantage  because  of  the  smaller  computational  effort  at  each  time  step. 

Explicit  direct  integration  schemes,  combined  with  mass  lumping  techniques,  have  been 
successfully  used  to  solve  dynamic  equilibrium  equations  for  the  single-phase  media.  The  main 
advantages  of  such  an  explicit  method  over  the  implicit  method  are  that  the  generalized  structural 
stiffness  matrix  does  not  have  to  be  formed  and  the  structural  equilibrium  equations  are  easily 
solved  by  just  inverting  the  diagonal  components  of  the  lumped  mass  matrix.  However,  the 
disadvantage  of  the  explicit  method  is  that  the  stable  time  step  size,  which  depends  on  the 
smallest  period  of  the  finite  element,  is  sometimes  unduly  small. 

This  section  presents  efficient  explicit  formulation  for  the  multi-phase  dynamic  analysis 
along  with  verification  problems  to  demonstrate  the  computational  algorithms  implemented  in 
the  Multi-Phase  Dynamic  Analysis  Program  (MPDAP). 
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Notation 


Note  that  positive  signs  have  been  used  for  elongation  and  tension.  A  comma  denotes 
differentiation  with  respect  to  the  subsequent  indices  and  the  superposed  dot  denotes 
difdferentiation  with  respect  to  time. 
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solid  phase  displacement 

apparent  fluid  displacement  relative  to  solid  skeleton 
total  stress 
effective  stress 
fluid  pressure 

solid  phase  displacement  vector  at  element  nodal  degrees  of  freedom 
apparent  fluid  displacement  vector  relative  to  skeleton  at  structural  nodal  degrees 
of  freedom 

applied  boundary  traction 

specified  boundary  pore  fluid  pressure 

elasto-plastic  stress-strain  matrix 

unit  vector  {1}^=  <  1  1  1000> 

porosity 

bulk  modulus  of  solid  grain 
bulk  modulus  of  soil/water  mixture 
bulk  modulus  of  elasto-plastic  skeleton 
bulk  mass  density  of  mixture 
fluid  mass  density 
dynamic  viscosity  of  pore  fluid 
linear  permeability  coefficient 
quadratic  permeability  coefficient 
Kronecker  delta 
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Explicit  Formulations 


The  fully  coupled  global  equilibrium  equations  for  the  multi-phase  medium  have  been 
derived  by  Blouin,  et.  al.  (1991)  and  are  given  by: 
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H  =  [N]^  [A^  dv 

=  E  I  K  ds 
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EE  =  {1}  {ir[5]  dv 

C  =  E|,^2  {^V[B]dv 
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The  inertial  force  terms  associated  mth  the  coupling  mass  matrix  and  in  Equation  2.4 
can  be  expressed  in  the  following  simple  form: 

(2.5) 

where: 


Employing  a  central  difference  scheme,  accelerations  and  velocities  at  step  n  can  be  expressed 
in  terms  of  displacement  increments  at  steps  n  and  n+1: 

and 

where  At  denotes  time  step. 


Now  substituting  Equations  2.6  and  2.7  into  2.5,  we  obtain: 


H  +^- M  •  A«  --J— H  •  Am 

2Ar  "  2Ar 


It  should  be  noted  that  for  practical  purposes,  the  mass  matrix  M  and  the  dissipation  matrix  H 
in  Equation  2.8  can  be  lumped  out  so  that  the  solution  can  be  performed  on  the  element  level 
without  the  factorization  of  the  structural  matrix. 
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Verification  of  Explicit  Algorithms 


Two  problems  are  presented  here  to  demonstrate  the  validity  of  explicit  algorithms 
implemented  in  MPDAP. 

The  first  problem  considers  a  12  inch  hollow  spherical  hole  in  an  infinite  elastic  medium. 
A  100  psi  step  load  is  suddenly  applied  to  the  internal  surface  of  the  spherical  hole  as  shown  in 
Figure  2.3.  Material  properties  and  time  steps  used  for  the  calculation  are  also  included  in 
Figure  2.3.  Figure  2.4  shows  the  radial  stress  profiles  at  5.5  msec.  The  MPDAP  calculation 
using  the  explicit  method  agrees  very  well  with  the  closed-form  solution. 

The  second  problem  considers  a  planar  compression  wave  propagating  through  saturated 
elastic  sand.  The  input  loading,  as  shown  in  Figure  2.5,  is  a  short  rise  time  triangular  pulse 
with  a  peak  stress  of  5,000  psi  and  a  positive  duration  of  10  msec.  Material  properties  and  time 
steps  used  for  calculations  are  also  included  in  Figure  2.5.  Two  MPDAP  calculations  were 
performed;  one  with  the  explicit  method  and  the  other  with  the  implicit  method. 

Figures  2.6  and  2.7  show  profiles  of  effective  vertical  stress  and  pore  pressure, 
respectively,  at  20  msec.  At  the  wavefront,  there  are  slight  differences  between  explicit  and 
implicit  methods.  But  behind  the  wavefront,  the  stress  profiles  are  almost  identical. 
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GIVEN:  Individual  Components  Properties  (Grain,  Pore-Fluid, 

and  Porous  Skeleton)  and  Specified  Mixture  Boundary 
Conditions/Load  Histories 
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Figure  2.1.  Material  element  model  outline. 
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Figure  2.2.  MEM  boundary  conditions. 


Time  Step  At  =  0.005  ms 

Young’s  Modulus  E  =  1 2,457  psi 
Poisson’s  Ratio  v  =  0.25 

Mass  Density  p  =  1 .88  x  1 0“*  Ib-s^/in^ 


Figure  2.3.  Verification  problem  1,  elastic  spherical 
wave  propagation  in  one-phase  medium. 
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Figure  2. 
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.  Radial  stress  profile  at  5.5  msec,  verification 
problem  1. 
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P  (psi) 


t(s) 


ASSUMED  MATERIAL  PROPERTIES 


Pore  Water 

0.29  X  10^  psi 

Solid  Grains 

Bulk  Modulus 

5.0  X  10*^  psi 

Specific  Gravity 

2.67 

Drained  Skeleton  Properties 

Bulk  Modulus 

3000  psi 

Constrained  Modulus 

6000  psi 

Poisson’s  Ratio 

0.20 

Porosity 

0.35 

Permeability 

1.0  in/s 

Explicit  Calculation 

At  =  0.8  X  10*  sec  and  C,  =  0.8 

Implicit  Calculation 

At  =  0.4  X  10^  sec,  iS  =  0.563 
and  7  =  1 

Figure  2.5.  Loading  time  history,  material  properties,  and 

time  steps  used  in  1-D  plane  strain  elastic  wave 
propagation  for  verification  problem  2. 
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SECTION  3 
MATERIAL  MODELS 


3.1  MODELING  APPROACH 

To  improve  our  understanding  of  the  mechanical  response  of  saturated,  porous  materials, 
we  performed  numerical  simulations  of  labcrutory  experiments  as  described  in  Sections  5  and 
7.  These  calculations  were  completed  using  the  finite  element  codes  MEM  and  MPDAP  which 
are  discussed  in  Section  2.  Central  to  the  success  of  these  calculations  are  the  numerical  models 
used  to  represent  the  materials.  These  models  are  described  in  this  section  with  the  application 
to  specific  material  given  in  later  sections  of  this  report. 

In  general,  our  approach  is  to  model  the  nonlinear  response  of  the  individual  material 
components;  that  is,  the  drained  stress-strain  response  of  the  porous  skeleton  and  the  pressure- 
volume  relationship  of  the  solid  grains,  pore  water,  and  pore  air.  The  behavior  of  these 
constituents  are  combined  using  fully  coupled  two-phase  theory  as  described  by  Kim,  Blouin, 
Chitty  and  Merkle  (1988)  and  Simons  (1989).  In  a  simple  decoupled  model,  the  soil  or  rock 
skeleton  acts  in  parallel  with  the  solid-fluid  mixture  such  that  the  net  compressibility  is  the  sum 
of  the  mixture  modulus  and  the  appropriate  skeleton  modulus.  The  shortcomings  of  a  decoupled 
model  are  demonstrated  by  Blouin  and  Kim  (1984).  A  fully  coupled  approach  incorporates  the 
mixture  model  with  the  additional  considerations  of  the  strain  in  the  skeleton  resulting  from  the 
compression  of  individual  grains  by  pore  pressure,  and  the  volume  change  in  the  solid-fluid 
mixture  due  to  effective  stresses  acting  on  the  solid  grains.  In  addition,  the  flow  of  pore  water 
within  the  porous  skeleton  (from  element  to  element  in  a  finite  element  mesh)  is  governed  by 
a  permeability  law.  This  sophisticated  formulation  provides  an  accurate  means  of  predicting 
undrained  material  response  for  many  porous  materials  up  to  very  high  pressures. 

This  section  describes  the  individual  components  of  the  material  model  in  detail.  Note 
that  not  all  of  these  models  were  used  in  this  study  but  are  included  for  completeness  in 
describing  the  models  available  in  the  MEM  and  MPDAP  codes.  These  models  include: 
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•  ARA  Three  Invariant  model  for  the  drained  response  of  the  porous  skeleton, 

•  nonlinear  solid  grain  model, 

•  pore  water  compressibility  model, 

•  partial  saturation  model  that  accounts  for  the  presence  of  pore  air,  and 

•  permeability  models. 


3.2  ARA  THREE  INVARIANT  SKELETON  MODEL 

In  simulating  the  behavior  of  saturated  materials,  it  is  particularly  important  to  accurately 
model  the  response  of  the  soil  skeleton.  Under  saturated  conditions,  the  material  strength  is 
strongly  influenced  by  the  pore  pressure  response  which  is,  in  turn,  controlled  by  the  volumetric 
deformation  characteristics  of  the  porous  skeleton.  That  is,  collapse  or  compressive  straining 
of  the  skeleton  transfers  more  load  to  the  fluid  in  the  pore  spaces  which  leads  to  reduced 
effective  stresses  and  shear  strength  in  the  material.  On  the  other  hand,  dilation  or  expansive 
strains  reduce  the  pore  pressures,  increase  the  effective  stresses,  and  generally  result  in  greater 
shear  strength.  Consequently,  the  skeleton  model  is  crucial  to  the  performance  of  any  saturated, 
two-phase  material  simulation. 

We  elected  to  use  a  relatively  advanced  skeleton  model  known  as  the  ARA  Three 
Invariant  (ARA  3-1)  model  which  was  originally  developed  by  ARA  under  sponsorship  of  the 
Air  Force  Office  of  Scientific  Research  and  described  by  Merkle  and  Dass  (1985)  and  Dass  and 
Merkle  (1986).  This  model  is  called  the  Three  Invariant  Model  because  the  controlling  surfaces 
in  stress  space  are  formulated  in  terms  of  three  independent  stress  invariants;  the  octahedral 
normal  stress,  octahedral  shear  stress,  and  Lode’s  angle.  Several  modifications  have  been  made 
to  permit  modeling  both  soil  and  rock  response  over  a  wider  range  of  stresses,  including  very 
high  stresses  associated  with  explosive  loadings. 

The  ARA  3-1  model  is  a  work  hardening,  elasto-plastic  model  with  two  independent  yield 
surfaces,  as  depicted  in  Figure  3.1.  The  compressive  yield  surface  is  an  ellipsoid  with  its  center 
at  the  origin  in  principal  stress  space.  It  is  associative,  only  strain  hardens,  and  the  strain 
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hardening  parameter  is  the  corresponding  plastic  work.  The  expansive  yield  surface  is  a 
hyperboloid  with  its  apex  on  the  hydrostatic  axis  in  principal  stress  space.  It  is  non-associative, 
strain  hardens,  and  the  strain  hardening  parameter  is  the  corresponding  plastic  work.  The 
parameters  for  the  ARA  3-1  model  are  summarized  in  Table  3.1. 

Fitting  and  evaluating  the  ARA  3-1  model  for  specific  materials  is  accomplished  using 
the  MEM  (Material  Element  Model)  code  outlined  in  Section  2.  One  portion  of  this  code 
contains  algorithms  for  fitting  the  model  while  another  part  is  used  to  exercise  the  model  under 
a  variety  of  boundary  conditions  on  a  single,  three  dimensional  finite  element.  The  fitting 
portion  of  the  MEM  code  outputs  the  ARA  3-1  model  parameters  in  the  format  described  in 
Table  3.2.  Other  codes  can  then  read  this  file  directly  and,  hence,  this  is  a  convenient  format 
for  presenting  the  model  parameters  for  a  given  material. 

The  model  is  fit  to  the  measured  response  of  the  material  skeleton.  Usually,  the  model 
is  fit  directly  to  laboratory  data  obtained  on  drained  test  specimens.  The  minimum  test  data 
required  to  fit  the  model  are: 

•  one  hydrostatic  (isotropic)  compression  test  with  at  least  two  unload/reload  cycles; 

•  three  triaxial  compression  tests  at  confining  pressures  in  the  range  of  interest  and 
with  unload  cycles  to  determine  Poisson’s  ratio; 

•  one  conventional  triaxial  extension  test  to  establish  the  material  strength  in  this 
mode. 

Other  test  data,  such  as  uniaxial  strain  tests,  undrained  triaxial  compression  tests,  etc.,  are  used 
to  evaluate  the  performance  of  the  model.  Section  5.1.1  of  this  report  presents  an  example  of 
fitting  the  ARA  3-1  model  to  test  data  for  a  porous  limestone. 
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3.2.1.  Notation 


sign  convention 

P. 

{e} 

{e'} 

{e/} 

Cv 

Tocl 

{o} 

^oct 

^oct 

*^oct,  max 
^ocl,  max,  TXC 

T 

J2’ 

J3’ 

(j) 

K 

Ko 

K, 

p 

fc 

fc' 

fc' 

fp 

fp 

fp' 

8p 


Compression  is  POSITIVE 
atmospheric  pressure 
total  strain  vector 
elastic  strain  vector 
plastic  compressive  strain  vector 
plastic  expansive  strain  vector 
volumetric  strain 
octahedral  shear  strain 
stress  vector 

octahedral  normal  stress  (mean  pressure) 

octahedral  shear  stress 

maximum  past  mean  pressure 

ultimate  octahedral  shear  stress  at  infinite  mean 

stress  in  triaxial  compression  mode 

tensile  strength  (negative) 

second  invariant  of  deviator  stress  tensor 

third  invariant  of  deviator  stress  tensor 

Lxide’s  angle 

bulk  modulus 

unload/reload  bulk  modulus  at  zero  pressure 
unload/reload  bulk  modulus  at  initial  unloading 
Poisson’s  ratio 
compressive  yield  criterion 
compressive  yield  stress  function 
compressive  yield  hardening  function 
expansive  yield  criterion 
expansive  yield  stress  function 
expansive  yield  hardening  function 
expansive  plastic  potential  function 
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We 

Wp 

pt 


V\ 

V2 

dXe 

dXp 

otf 


compressive  plastic  work 

expansive  plastic  work 

expansive  plastic  work  at  peak  shear  stress 

normalized  expansive  plastic  work 

ratio  between  the  triaxial  compression  and  extension 

strengths  at  a  given  mean  stress 

expansive  yield  hardening  parameter 

expansive  plastic  potential  parameter 

compressive  yield  proportionality  constant 

expansive  yield  proportionality  constant 

fraction  of  ultimate  in  triaxial  compression  mode 


3.2.2.  Total  Strain  Formulation 

The  ARA  Three  Invariant  model  computes  the  total  strain  in  parts  consisting  of  the 
elastic,  compressi'-e  plastic,  and  expansive  plastic  strain  components; 

{de}  =  {de^}  *  {dd>J  +  {deP^}  (3.1) 


Each  strain  component  is  calculated  independently  within  the  model  as  described  in  the  following 
subsections. 

3.2.3.  Elastic  Response 

At  stress  states  inside  the  yield  surfaces,  the  skeleton  response  is  treated  as  nonlinear 
elastic  and  governed  by  the  previous  maximum  peak  stress.  Two  options  are  available  for 
modeling  the  elastic  response  within  the  framework  of  the  three  invariant  skeleton  model:  the 
ARA  elastic  model  and  the  i.  '.e  and  Nelson  elastic  model.  In  both  options,  Poisson’s  ratio  is 
assumed  to  remain  constant. 
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ARA  elastic  model.  During  virgin  loading,  the  elastic  bulk  modulus  is  given  by: 


K  = 


3  (1  -  20 


>  K, 


(3.2) 


where  K^,  and  n  are  material  constants  obtained  in  the  parameter  fitting.  K;  represents  the  initial 
bulk  modulus  at  low  pressures  and  is  necessary  for  modeling  the  behavior  of  rock-type  materials 
that  have  a  definite  initial  elastic  behavior.  In  uncemented  soils,  can  be  taken  as  a  very  small 
value.  The  initial  bulk  modulus  is  also  used  to  determine  the  initial  position  of  the  compressive 
plastic  yield  function  (cap)  by  defining  the  initial  elastic  range. 


During  unloading  or  reloading,  the  skeleton  modulus  is  described  by  one  of  two  segments 
as  depicted  in  Figure  3.2.  Between  the  previous  peak  mean  stress,  and  the  transition 

into  the  nonlinear  segment  at  (7«t,b>  the  elastic  bulk  modulus  is  constant  and  is  given  by; 


K  =  K, 


3  (1  -  2u) 


^oci.nutx 


n 


The  transition  into  the  nonlinear  segment  occurs  at: 

^oct.b  ^  ^or/.nux 


(3.3) 


(3.4) 


where  X  is  a  model  parameter.  At  mean  stresses  less  than  the  nonlinear  bulk  modulus  is 
given  by; 


(3.5) 


Referring  to  Figure  3.2,  the  model  parameters  7  and  /3  are  given  by: 


where  is  the  bulk  modulus  at  zero  pressure  and 


(3.6) 
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(3.7) 


K 


where  K.  is  the  bulk  modulus  at  one  quarter  of  the  transition  pressure  5. 


While  this  formulation  allows  for  relatively  accurate  curve  fitting  of  observed  soil 
response,  the  model  has  three  disadvantages: 

a)  for  certain  closed-loop  stress/strain  paths,  the  model  may  violate  the  energy 
conservation  principle; 

b)  unloading  at  low  pressures  could  potentially  generate  expansive  volumetric 
strains;  and 

c)  at  the  transition  pressure,  the  modulus  is  not  continuous. 


Lade  and  Nelson  elastic  model.  The  second  elastic  model  option  is  based  on  a 
relationship  derived  by  Lade  and  Nelson  (1987).  This  formulation  is  continuous  and  was  derived 
from  the  energy  conservation  principle.  Lade  and  Nelson’s  model  can  be  expressed  as: 


-  .. 

2 

n 

7 

KP. 

3  <^oc, 

^  6  (1  +  J^)  -^2 

3"'*  (1  -  2u) 

Pa 

1  -  2p  p2 

0 

>  K.. 


(3.8) 


where  the  parameters  n,  and  Kj  are  the  same  as  used  in  Equation  3.2.  Since  this  model  is 
fit  strictly  using  the  slope  of  an  initial  unload  curve,  it  can  be  difficult  to  closely  match  the 
observed  characteristic  of  an  unload  cycle. 

Fitting.  Poisson’s  ratio  (i')  for  a  given  material  can  be  determined  in  a  number  of  ways 
using  unload/ reload  data  which  represents  the  elastic  response  of  the  skeleton.  Lade  and  Nelson 
(1987)  recommend  obtaining  Poisson’s  ratio  directly  from  strain  measurements  in  triaxial 
compression  unload/reload  cycles  right  after  stress  reversal  at  hydrostatic  conditions  where: 
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where 

e,  =  axial  strain 

e,  =  radial  strain 

ev  =  volume  strain 

(3.9) 


In  addition,  since  the  elastic  response  is  completely  defined  by  any  two  independent 
elastic  parameters,  Poisson’s  ratio  can  be  obtained  from  the  bulk  modulus,  (K),  measured  in 
a  hydrostatic  compression  unload,  and  any  other  elastic  modulus.  For  example,  a  triaxial 
compression  unload  yields  the  shear  modulus  (G),  an  unconfmed  compression  unload  gives  the 
Young’s  modulus  (E),  and  a  uniaxial  strain  unload  produces  the  constrained  modulus  (M).  Any 
one  of  these  parameters  can  be  used  with  the  bulk  modulus  to  obtain  Poisson’s  ratio; 


3K  -  2G 
2  (3K  +  G) 


(3.10) 


3K  -  E 
6K 


(3.11) 


3K  -  M 
3K  +  M 


(3.12) 


To  obtain  the  elastic  model  parameters  and  n.  Equation  3.2  is  rewritten  in  the  form: 


log 


3K  (1  -  2j^) 
Pa 


=  log  Kr  *  O  log 


(3.13) 


Values  of  K  and  (Toc,,  from  the  initial  unloading  response  at  various  pressures  in  the  hydrostatic 
compression  test,  are  then  plotted  as  log  (3K  (l-2v)/PJ  versus  log  (  A  least  squares 

linear  regression  is  then  applied  in  log-log  space.  The  parameter  n  is  the  slope  of  this  line  while 
Ku,  is  the  intercept  where  (<roct/PJ  is  1.0.  The  parameters  X,  y,  and  for  the  ARA  elastic 
unload  model  are  determined  from  a  single  unload/reload  cycle  in  the  hydrostatic  compression 
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test  as  depicted  in  Figure  3.2.  The  parameters  are  computed  using  Equations  3.4,  3.6,  and  3.7. 


3.2.4  Compressive  Plastic  Response 

The  compressive  yield  surface,  or  cap,  is  an  ellipsoid  with  its  center  at  the  origin  in 
principal  stress  space  as  shown  in  Figure  3.1.  The  ratio  of  the  major  axis  to  the  minor  axis  is 
defined  by  the  constant  r,  as  shown  in  Figure  3.1,  and  is  determined  in  the  parameter  fitting 
process.  The  compressive  yield  criterion  is  given  by 

fc  =//  ■//  =  0 

where  the  compressive  yield  stress  function,  f/,  is  given  by 

f  =  3(ff  2  +  r  (3.15) 

and  the  compressive  yield  hardening  function,  f/',  is  given  by 

(3.16) 


r 


In  the  model,  f/'  is  defined  in  segments  with  the  parameters  Cj  and  p;  determined  for  each 
segment  in  the  fitting  process.  W,  is  the  compressive  plastic  work  given  by 

The  compressive  flow  rule,  which  is  associative,  is  given  by 


{*/)  • 


do 


(3.18) 


where  d\.  is  the  compressive  yield  proportionality  constant. 
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Fitting.  The  elliptical  cap  ratio,  r,  can  be  determined  directly  from  special  tests 
following  a  hydro-triax-hydro  type  stress  path  (Applied  Research  Associates,  Inc.,  1991). 
However,  this  special  test  data  is  seldom  available  and  the  results  are  still  subject  to  considerable 
interpretation.  Hence,  the  cap  ratio  is  generally  selected  to  give  reasonable  model  behavior 
based  on  experience  and  engineering  Judgement. 


Other  parameters  which  define  the  compressive  plastic  response  are  derived  directly  from 
the  hydrostatic  compression  test  data.  At  a  given  pressure  in  the  hydrostatic  data,  the  elastic 
component  of  the  volume  strain  is  computed  from  the  elastic  model  and  parameters  described 
in  Section  3.2.3.  The  compressive  plastic  strain  is  obtained  by  subtracting  the  elastic  component 
from  the  total  measured  volume  strain  and  the  compressive  plastic  work  is  then  computed 
according  to  Equation  3.17.  Then,  for  hydrostatic  compression,  the  octahedral  shear  stress  is 
zero  and  Equations  3.14,  3.15,  and  3.16  reduce  to; 


W 


C 


*  C; 


3 


(3.19) 


Taking  the  logarithm  of  both  sides  yields: 


log 


-  .. 

_ 

2 

w 

a  , 

c 

Pa 

=  log  C.  +  p.  log 

3 

oct 

'K 

(3.20) 


Plotting  log  (Wj/PJ  versus  log  [3  (ffoc/PJ^l  allows  us  to  determine  the  parameters  C;  and  p;  as 
the  intercepts  and  slopes,  respectively,  of  straight  line  fits  to  the  data.  Since  the  shape  of  this 
data  is  usually  curved  in  log-log  space,  multiple  straight  line  segments  are  defined  by  fitting  q 
and  Pi  over  specific  ranges. 
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3.2.5  Expansive  Plastic  Response 


Yield  surface.  The  expansive  yield  surface,  or  shear  yield  surface,  is  a  hyperboloid  with 
its  apex  on  the  hydrostatic  axis  in  principal  stress  space  as  shown  in  Figure  3.1.  The  expansive 
yield  surface  has  the  same  shape  as,  and  approaches  asymptotically  to,  the  ultimate  failure 
surface  as  shown  in  the  same  figure.  The  shape  of  the  expansive  yield  surface  in  the  -r  plane, 
perpendicular  to  the  hydrostatic  axis  (see  Figure  3.1),  is  a  triple  ellipse  in  polar  coordinates. 


The  expansive  yield  criterion  is  given  by: 

4  =/;  -//  =  0 


where  the  expansive  yield  stress  function,  fp'  is  given  by 


' 

1  -  cos  3a) 

\  E  j 

(3.21) 


(3.22) 


where  E  and  m  are  parameters  determined  in  the  fitting  process. 


E  is  computed  as: 


E  rP  ^  I 


(3.23) 


where  'ir  is  the  ratio  between  the  triaxial  compression  and  extension  strengths  at  a  given  mean 
stress  as  indicated  in  Figure  3.1.  It  should  be  noted  that  ^  must  be  less  than  9/7  or  the  ultimate 
failure  surface  will  become  convex  in  the  ir  plane  thereby  violating  uniqueness.  Lode’s  angle, 
w,  shown  in  Figure  3.1,  is  defined  in  the  t  plane  and  varies  from  0°  to  120°.  An  cj  of  120° 
corresponds  to  the  triaxial  compression  mode  (coi'  3w  =  1.0)  and  a  value  of  60°  corresponds 
to  the  triaxial  extension  mode.  Lode’s  angle  is  related  to  the  second  and  third  deviatoric  stress 
invariants  by: 
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cos  3aj  = 


1.5 


(3.24) 


The  expansive  yield  surface  hardens  isotopically  with  the  yield  hardening  function: 

/;' -1,  <3.25) 


The  coefficient  rji  is  a  failure  parameter  representing  the  maximum  value  of  fp'  and  q  is  an 
expansive  hardening  parameter  determined  in  the  parameter  fitting  process,  is  the  normalized 
expansive  plastic  work  given  by 


W 

P 

f 

p.peak 


(3.26) 


where  the  expansive  plastic  work,  Wp,  is  given  as 

»  I  {</(/)  (3.27) 

and  Wp  pejk  is  the  peak  plastic  expansive  work  at  the  ultimate  failure  envelope  for  a  given  mean 
stress. 


In  the  model,  two  options  exist  for  defining  Wp  p^^ik.  The  first  option  uses  a  multilinear 
function  of  mean  stress: 


W 


p,peak 


=  o,  <r„,  +  b, 


(3.28) 


The  coefficients  and  b;  are  determined  in  the  parameter  fitting  process  for  specific  ranges  of 
mean  stress.  This  option  was  added  to  the  model  to  permit  more  accurate  modeling  of  rock. 
The  second  option  for  defining  Wp  p^.^  uses  an  exponential  form; 
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w 


p,ptak 


-P  Pa 


(3.29) 


where  p  and  /  are  model  parameters. 


Triaxial  compression  strength.  In  evaluating  a  model  fit,  it  is  useful  to  plot  the  strength 
envelope  for  triaxial  compression  stress  paths.  On  the  strength  envelope,  the  plastic  work  equals 
Wppj^  such  that  |p  =  1.0.  Then,  Equation  3.25  becomes 

//  =  77,  (1.0)  =  77,  >'3.30) 


In  triaxial  compression  mode,  cos  3u  equals  1.0  and  Equation  3.22  reduces  to: 


(3.31) 


Substitution  of  Equations  3.30  and  3.31  into  Equation  3.21  yields  an  equation  for  the  triaxial 
compression  strength  envelope  in  the  octahedral  stress  space: 


oct  _ 


-T 


1  +  m 


<^<.7  -T 


1  - 


(3.32) 


As  the  mean  stress  goes  to  infinity,  the  triaxial  compression  strength  envelope  asymptotically 
approaches  the  ultimate  shear  strength  given  by 


^ocl,ma,TXC  _  1 

K  m 


1  - 


(3.33) 
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Expansive  flow  rule.  The  expansive  flow  rule,  which  governs  the  nonassociative 
expansive  plastic  flow,  is  given  by 


(3.34) 


where  dXp  is  the  expansive  yield  proportionality  constant.  The  incremental  expansive  plastic 
strain  vector  is  perpendicular  to  the  expansive  plastic  potential  function  defined  by 


S,  = 


V2 


(1  -  g  cos  3o)) 


-T 


(3.35) 


1  +  m 


-T 


The  plastic  potential  parameter  tj,  can  be  defined  in  two  ways  within  the  model.  The  first 
option  defines  the  parameter  in  terms  of  a  multi-segment  function  of  the  expansive  yield 
hardening: 

(3-36) 

where  tj  and  S;  are  determined  for  specific  ranges  of  fp*  in  the  fitting  process.  This  option  was 
added  to  the  model  to  improve  the  capability  for  modeling  rock.  The  second  option  employs 
a  nonlinear  function  of  the  form: 


1\2  =  t  +  R 


_oa  +  s 

P  ''P 


(3.37) 


where  the  parameters  t,  R,  and  s  are  constant  at  all  ranges.  This  option  produces  a  continuous 
function  for  tj2. 


Potential  function  at  high  pressure.  In  general,  the  ARA  3-1  model  uses  nonassociative 
flow  for  computing  plastic  strain  increments  related  to  the  expansive  or  shear  yield  surface. 
That  is,  the  expansive  plastic  strain  vector,  {ep’’},  is  computed  perpendicular  to  the  plastic 
potential  function  given  in  Equation  3.35.  However,  in  the  special  case  where  172  =  fp’, 
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Equation  3.35  can  be  shown  to  reduce  to: 


8p 


(3.38) 


Thus,  for  the  special  case  when  r]2  =  fp',  the  potential  function  is  some  fraction  of  the  expansive 
yield  surface  and  plastic  flow  becomes  associative. 


At  very  high  pressures,  where  the  expansive  yield  surface  is  approaching  the  ultimate 
strength  asymptote,  the  material  should  behave  like  a  fluid  where  shear  stresses  do  not  generate 
plastic  volume  strains.  To  model  this  high  pressure  behavior,  the  plastic  potential  function  is 
modified  in  the  higher  pressure  regime  to  yield  associative  flow,  as  depicted  in  Figure  3.3. 
First,  we  define  otf  as  the  fraction  of  the  ultimate  octahedral  shear  stress  at  infinite  mean 
pressures: 


T 


OCt,<Xf 


oct.tD*x.TXC 


(3.39) 


Then,  using  the  strength  envelope,  the  corresponding  mean  stress,  cr„„  can  be  found  for  any 
given  fraction  oif. 


As  indicated  in  Figure  3.3,  the  potential  function  is  defined  in  three  pressure  regimes 
delineated  by  Of  equal  to  0.5  and  0.8.  In  the  low  pressure  regime  (a*,  <  (Tocuof  =  5).  the  potential 
function  parameter  tjj  is  determined  using  the  values  of  s,  t,  and  R  derived  from  the  model  fit. 
At  very  high  pressures  (a^ct  >  (Toct,  of  =  o.»)>  the  tjj  parameters  are  arbitrarily  set  to  produce 
associate  flow.  Finally,  in  the  transition  zone  {Ogcu  os  =  o.s  <  ‘^oct  <  =  o.s)  the  s,  t,  and  R 

parameters  are  linearly  interpolated. 

Fitting.  The  first  step  in  determining  the  expansive  plastic  response  parameters  is  to 
define  the  ultimate  failure  surface  and  the  corresponding  peak  expansive  plastic  work  as 
functions  of  mean  stress.  In  the  conventional  manner,  the  failure  envelope  is  fit  to  the  triaxial 
compression  test  data  at  the  peak  shear  stress  or  the  shear  stress  at  15  %  axial  strain  if  no  peak 
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is  observed.  The  tensile  strength  T  can  be  measured  directly  in  an  unconfmed  tensile  test  (for 
rocks)  or  inferred  from  the  shape  of  the  failure  surface.  For  uncemented  soils,  T  is  generally 
very  small.  In  the  triaxial  compression  mode,  the  failure  surface  is  defined  by  Equation  3.32 
which  can  be  rearranged  to  the  form: 


' 

-  T  _  1  “  E  ^  ^ 

^  ■  E 

-T 

7i 

Pa 

(3.40) 


By  plotting  the  failure  stress  points  from  each  triaxial  compression  test  in  terms  of  -  T)/r<^, 
versus  (ffo^,  -  T)/P„  a  straight  line  fit  will  yield  an  intercept  of  (1  -  E)/rj,  and  a  slope  of 
[m  (l-E)/rj,].  Then,  the  parameter  m  is  obtained  simply  as  the  ratio  of  the  slope  to  the  intercept 
of  this  line.  The  parameter  E  is  obtained  directly  from  the  ratio  of  triaxial  compression  and 
triaxial  extension  strength  (’^^)  as  defined  in  Equation  3.23.  Substitution  of  E  into  the  intercept 
value  yields  the  parameter  jj,.  This  completes  the  definition  of  the  shear  failure  surface. 


The  next  step  is  to  determine  the  peak  plastic  expansive  work  as  expressed  in  Equations 
3.28  or  3.29.  The  plastic  expansive  strain  is  first  computed  by  subtracting  the  elastic  and 
compressive  plastic  strain  components  from  the  total  strains  in  the  triaxial  compression  tests. 
The  plastic  expansive  work  at  the  peak  stress,  Wppjjk,  is  computed  by  integration  according  to 
Equation  3.27.  If  Wp  is  to  be  modeled  with  Equation  3.28,  the  values  are  plotted  versus  the 
mean  stress  and  fit  with  straight  line  segments.  On  the  other  hand,  if  Wpp*^  is  to  be  modeled 
using  Equation  3.29,  we  first  rewrite  the  relationship  as: 

=  log  p  +  I  log 


OCl 

T. 


(3.41) 


A  least  squares,  linear  fit  to  log  (Wp  p^^^/PJ  versus  log  will  then  yield  1  and  p  from  the 

slope  and  intercept,  respectively. 

The  next  step  is  to  define  the  expansive  yield  hardening  function  of  Equation  3.25.  The 
coefficient  ijj  has  been  obtained  previously  from  the  fit  to  Equation  3.40  and  the  peak  plastic 
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expansive  work  is  defined  by  Equation  3.28  or  3.29.  Thus,  only  the  exponent  q  must  be 
obtained  in  order  to  completely  define  the  hardening  function  fp'.  Solving  Equation  3.25  for  q 
gives: 


(J  = 


In 

/" 

jp 

7, 

(1 

plot  of  fp'Vrj,  as 

(3.42) 


triaxial  test  data  at  each  confining  pressure.  We  then  match  the  hardening  function  for  each 
triaxial  test  at  a  value  of  fp'Vjji  equal  to  0.80  according  to  Equation  3.42.  The  q  values  obtained 
in  this  manner  are  then  averaged  to  produce  the  parameter  q  for  the  model  fit. 


The  final  step  remaining  is  to  specify  the  expansive  flow  rule  and  expansive  plastic 
potential  function  as  defined  in  Equations  3.34  and  3.35  and  illustrated  in  Figure  3.1.  The 
expansive  plastic  strain  increment  is  perpendicular  to  the  plastic  potential  function  as  defined  in 
Equation  3.35.  Below  50%  of  the  ultimate  shear  strength,  the  potential  function  has  a  form 
similar  to  the  expansive  yield  function  given  in  Equation  3.22,  but  in  general  the  two  functions 
are  not  equal  and  the  expansive  plastic  strain  is  not  normal  to  the  yield  surface;  i.e.,  a 
nonassociative  flow  rule  is  used.  Recall  from  the  previous  subsection  that  the  potential  function 
transitions  to  give  associative  flow  at  very  high  pressures  as  depicted  in  Figure  3.3.  In  order 
to  define  the  plastic  potential  function,  we  need  only  to  determine  the  parameter  tjj,  since  all 
other  parameters  have  been  previously  derived. 

The  physical  meaning  of  parameter  rjj  is  depicted  in  Figure  3.4  at  the  intersection  of  the 
potential  surface  and  the  compressive  yield  surface.  The  expansive  plastic  strain  vector  (dip’’), 
compressive  plastic  strain  vector  (di/)  and  total  plastic  strain  vector  (di^)  are  depicted  at  the 
point  of  intersection.  From  normality,  the  ratio  of  the  shear  component  of  the  expansive  plastic 
strain  vector  to  the  volumetric  component  equals  the  ratio  of  the  partial  derivative  of  the 
expansive  potential  function  with  respect  to  the  shear  stress  to  the  partial  derivative  of  the 
potential  function  with  respect  to  the  mean  stress  as  given  by: 
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(3.43) 


3  ^ 


Substituting  Equation  3.35  for  in  Equation  3.43  and  solving  for  tjj  gives: 


^2 


_1 

(i  -  -] 

a^,-T 

3 

^yoc.p 

\  Ej 

Pa 

(3.44) 


In  the  triaxial  compression  test; 


oct.p 


4 


(3.45) 


and 

*  2*V 

where  e,  and  e,  denote  the  axial  and  radial  strains,  respectively.  Equation  3.44  shows  that  the 
parameter  1J2  is  linearly  proportional  to  the  inverse  slope  of  the  plastic  potential  surface. 

Hence,  the  plastic  potential  function  parameter  is  computed  at  points  along  the  triaxial 
compression  response  using  Equation  3.44.  Then,  if  Equation  3.36  is  employed,  multi-segment 
linear  fits  are  determined  for  tjj  as  a  function  of  fp'.  If  Equation  3.37  is  used  to  define  tjj  ,  a 
regression  analysis  is  used  to  determine  t,  R,  and  s  for  a  single  function  defining  tjj. 

This  completes  specification  of  the  ARA  Three  Invariant  model  and  derivation  of  all 
model  parameters  from  the  test  data.  A  summary  of  the  model  parameters  is  given  in  Table  3.1. 

3.2.6  Tensile  Failure  of  Skeleton 

It  is  possible,  particularly  in  materials  around  an  explosive  detonation,  for  the  skeleton 
to  fail  in  tension.  In  a  soil  or  rock  mass,  the  material  is  literally  pulled  apart  such  that 
individual  grains  are  physically  separated  or  pieces  of  unconnected  material  may  form.  In 
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modeling  rocks  or  concrete,  this  is  a  spall  condition  and  occurs  when  the  tensile  strength  is 
exceeded.  After  failing  in  tension,  the  skeleton  has  no  stiffness  but  can  collapse  back  together 
to  again  form  a  contiguous  mass.  While  the  exact  nature  of  this  behavior  is  not  easily  described 
for  the  general  case,  particularly  the  characteristics  of  a  reformed  skeleton,  the  material  model 
must  be  able  to  track  the  appropriate  skeleton  behavior  under  these  conditions. 

In  the  ARA  3-1  skeleton  model,  the  tensile  strength  is  given  by  the  parameter  T  which 
is  the  apex  of  the  expansive  plastic  strength  envelope.  When  tensile  stresses  exceed  the  tensile 
strength,  the  volume  strain  at  which  the  tensile  failure  occurred  is  stored.  The  skeleton  stiffness 
is  set  to  zero  and  the  three  principal  stresses  are  set  equal  to  the  tensile  strength.  The  volume 
strains,  which  are  initially  expansive,  are  then  tracked  with  respect  to  the  volume  strain  at  tensile 
failure. 


At  some  point  the  volume  strains  may  become  compressive  and  again  return  to  the  total 
volume  strain  at  which  the  tensile  failure  occurred.  At  this  point,  the  material  is  considered  to 
be  reformed  and  possess  the  same  stiffness  and  properties  as  the  soil  had  when  the  tensile  failure 
first  occurred.  In  reality,  we  would  expect  the  reformed  skeleton  to  be  significantly  different 
from  when  the  tensile  failure  occurred.  While  this  model  does  not  consider  this  altered 
behavior,  it  does  ensure  that  the  skeleton  will  begin  to  support  compressive  loads  upon 
reforming. 


3.3  GRAW  COMPRESSIBILITY  MODEL 

To  model  the  nonlinear  response  of  the  solid  grains  to  both  the  applied  pore  pressure  and 
effective  stress,  analytic  expressions  for  the  deformation  of  solids  at  high  pressure  are  employed. 
High  pressure  data  for  many  rocks  and  minerals  show  a  linear  relationship  between  loading  wave 
velocity  and  particle  velocity  (e.g.  Allen,  1967).  The  loading  wave  velocity  can  be  expressed 
as: 
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where: 


Cl  =  loading  wave  velocity 

Co  =  the  initial  wave  velocity  at  relatively  low  pressure 
Vp  =  peak  particle  velocity 

S  =  experimentally  determined  constant  relating  Cl  to  Vp  (generally 
equal  to  about  1.5  for  most  dense  rocks  and  minerals). 


Conservation  of  mass  and  momentum  on  either  side  of  the  wavefront  yield  the  familiar 
relationships: 


a  =  p  c 

p  i 


(3.48) 


where: 


<Tp  =  peak  axial  stress 
Po  =  initial  material  density 
M  =  constrained  secant  modulus  = 

€p  =  peak  axial  strain  corresponding  to  the  peak  stress 


Substitution  of  Equation  3.47  into  3.48  gives: 

0  =  p  c  V  +  p„Sv^ 

p  o  p  '^o  p 


(3.50) 


and  solving  for  peak  particle  velocity  as  a  function  of  peak  stress  yields 

V  =  (3-51) 

2p,S 


where: 


Substitution  of  Equation  3.47,  3.51,  and  3.52  into  Equation  3.49  gives: 


M  =  F(ap  =  +  cjia)  +  (3.53) 

The  tangent  constrained  modulus,  M„  used  in  the  numerical  model  is  defined  as  the  slope 
of  the  stress  strain  curve  by: 


da 


(3.54) 


From  Equation  3.53  and  the  definition  of  constrained  modulus,  M: 


e 


P 


(3.55) 


Differentiating  Equation  3.55  with  respect  to  a^  and  inverting  gives  the  tangent  constrained 
modulus  as 


M.  = 


-  (^pF'io^) 


(3.56) 


Differentiating  Equations  3.52  and  3.53  with  respect  to  a^  yields: 


F'{o;)  =  cj'ia^)  * 

^Po 


(3.57) 


and 


= 


(P.V  ^  4p>p>'^ 


(3.58) 


Hence,  Equations  3.51  through  3.58  can  be  used  to  define  the  high  pressure  constrained  stress 
strain  and  modulus  relationships  for  the  solid  grains. 


For  two  phase,  coupled  calculations,  the  volumetric  relationships  for  the  solid  grains 
should  be  specified  in  terms  of  the  bulk  modulus.  Kg,  rather  than  in  terms  of  the  constrained 
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modulus.  At  high  pressures,  the  shear  strength  of  the  grain  materials  becomes  insignificant 
compared  to  the  applied  stress  and  the  materials  tend  to  behave  like  fluids.  At  these  pressures, 
the  tangent  bulk  modulus  equals  the  tangent  constrained  modulus  with  Poisson’s  ratio  equal  to 
0.5.  Beneath  some  threshold  pressure,  pb,  Poisson’s  ratio  begins  to  decrease  from  0.5  at  p,,  to 
an  initial  value  of  Poisson’s  ratio,  at  a  low  value  of  mean  stress.  We  have  used  a  simple 
relationship  to  approximate  the  influence  of  mean  stress  on  Poisson’s  ratio  for  the  solid  grains: 

=  8(PW,  (3.59) 


The  ratio  of  the  bulk  modulus  to  the  tangent  constrained  modulus,  g(p)  at  pressures  less  than  pt 
is  given  by: 


2  (1  -  ^  ((  ^  O 

3  (1  -  O  Pb  3(1  -  v) 


(3.60) 


For  pressures  greater  than  p^; 

g{p)  =  1  (3.61) 

Poisson’s  ratio  can  be  computed  as  a  function  of  the  modulus  ratio  at  a  given  pressure  as: 

^  3^(P)  -  1  (3.62) 

I  'igip) 


3.4  PORE  WATER  COMPRESSIBILITY 

The  model  for  the  nonlinear,  elastic  compressibility  of  the  pore  water  is  derived  from  an 
equation  of  state  reported  by  Ahrens  (1988)  and  attributed  to  Bakanova,  et.  al.  (1976).  This 
equation  relates  the  shock  velocity  in  water  to  the  peak  particle  velocity.  In  the  lower  pressure 
regime,  a  quadratic  relation  is  used  while  a  linear  relation  is  used  in  the  higher  pressure  regime. 
The  transition  point  between  the  two  regimes  is  defined  in  terms  of  a  peak  particle  velocity  at 
the  transition,  Vp,.  Bakanova’ s  equations  can  be  expressed  as: 
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(3.63) 


V  <  V  .  ■ 

Vp  —  »pj  . 


Vp  >  Vp.  : 


c  =  Cj  + 


where:  c 

Vp 

Ci>  Si,  S2 
^2.  S3 


shock  propagation  velocity  in  the  fluid 
peak  fluid  particle  velocity 
constants  used  to  fit  data  below  the  transition 
constants  used  to  fit  data  above  the  transition 


(3.64) 


Equation  3.64  can  also  be  expressed  in  terms  of  the  shock  velocity  at  the  transition  point,  c,. 
Substituting  Vp,  into  Equation  3.64  yields: 


Substituting  3.65  into  3.64  produces  this  expression  for  the  shock  velocity  above  the  transition: 


Vp  >  Vp,: 


c  =  c,  +  53  (V  -  V  ) 


(3.66) 


where; 


c,  =  shock  velocity  at  the  transition 

Vp,  =  peak  particle  velocity  at  the  transition  (model  constant) 


At  the  transition  point,  the  shock  velocity  from  Equations  3.63  and  3.66  should  be  equal  to 
preserve  continuity.  Setting  Equations  3.63  and  3.66  equal  at  Vp  =  Vp,  gives: 


c  2 


(3.67) 


thereby  defining  c,  in  terms  of  the  model  constants.  Equations  3.63,  3.66,  and  3.67  (with  the 
constants  c,,  Sj,  S2,  and  S3)  define  the  shock  velocity  as  a  function  of  peak  particle  velocity. 

To  derive  a  bulk  modulus  for  water  as  a  function  of  pressure,  we  first  need  an  expression 
for  the  peak  particle  velocity  as  a  function  of  pressure.  Conservation  of  mass  and  momentum 
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on  either  side  of  the  wavefront  yields  the  familiar  relationship  from  shock  physics: 


ir„  =  pcv^ 

P  p 


(3.68) 


where: 


Tp  =  pore  fluid  pressure 
Po  =  mass  density  of  fluid 


Substitution  of  Equation  3.63  into  3.68  yields  an  expression  for  the  transition  fluid  pressure  (Xpi): 


(3.69) 


For  water,  the  transition  pressure  is  greater  than  30,000  MPa. 


Below  the  transition  pressure,  substitution  of  Equation  3.63  into  3.68  will  give: 


+  ^  V  -  =  0 


(3.70) 


This  cubic  equation  can  be  solved  to  yield  an  expression  for  Vp  as  a  function  of  fluid  pressure 
below  the  transition  pressure  iTp^ 


Vp  =  m  cos 


1  f  3/3  1  .  4ir1  _ 

_  cos  -C—  +  _—  - 

3  am  3  J  35^ 


(3.71) 


where: 


a  =  £l  -  I 
“  "  5^  3  52 


(3.72) 


p  5,  3  5,  5,  27  5, 


(3.73) 


~oi 

m  =  2  — 


(3.74) 
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Above  the  transition  pressure,  substitution  of  Equation  3.66  into  3.68  yields  a  quadratic 
equation: 


(3.75) 


Solving  this  equation  for  Vp  as  a  function  of  fluid  pressure  gives  Vp  for  pressures  above  the 
transition  pressure 


’ 

2 

253 

253 

P/3_ 

(3.76) 


The  elastic  bulk  modulus  of  water  (K,„)  is  defined  as: 

d-Kjdv^ 

P  P 


(3.77) 


where  Cy  is  the  volume  strain  corresponding  to  the  pressure  Tp.  Taking  the  derivative  of 
Equation  3.68: 


The  volume  strain  is  given  by: 


€ 


P 

c 


and  taking  the  derivative  yields: 


^  ^  ^  ~  V' 
dv^ 


(3.78) 


(3.79) 


(3.80) 


Substitution  of  Equations  3.78  and  3.80  into  3.77  gives  an  expression  for  the  bulk  modulus  in 
terms  of  the  shock  and  peak  particle  velocities: 
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(3.81) 


= 


c  -  v^c 


J 


The  derivatives  of  the  shock  velocity  with  respect  to  the  peak  particle  velocity  are  given  by; 

c'  =5,  (3.82) 


For  Tp  <  Tp,; 


For  Tp  >  Tp,; 


(3.83) 


This  completes  the  derivation  of  the  water  compressibility  model.  To  summarize  the  use 
of  the  model: 

1.  begin  by  computing  the  transition  shock  velocity  and  pressure  using  Equations 
3.67  and  3.69, 

2.  compute  the  peak  particle  velocity  using  either  Equation  3.71  below  the  transition 
pressure  or  Equation  3.76  above  the  transition  pressure, 

3.  compute  the  corresponding  shock  propagation  velocity  from  Equation  3.63  or 
3,66  and  the  derivatives  from  Equation  3.82  or  3.83; 

substitute  the  results  into  Equation  3.81  to  get  the  bulk  modulus. 


The  material  constant  values  for  this  model  are  given  in  Table  3.3  for  fresh  water  and  sea  water. 
The  fresh  water  values  are  from  Bakanova,  et.  al.  (1976)  as  reported  by  Ahrens  (1988). 
Parameters  for  sea  water  were  fit  to  compressibility  data  described  by  Kim,  et.  al.  (1986)  and 
attributed  to  Britt  (1985).  The  pressure- volume  response  of  fresh  water  and  sea  water  as 
predicted  by  this  model  are  given  in  two  pressure  ranges  in  Figures  3.5a  and  3.5b. 


Distention  and  Rejoining.  Just  as  it  is  possible  for  the  skeleton  to  separate  into  pieces 
in  a  tensile  failure,  gas-filled  voids  develop  in  the  pore  fluid  if  subjected  to  tensile  (vacuum) 
pressures.  These  voids  will  contain  either  water  vapor  or  air  drawn  in  from  the  atmosphere  and 
may  exist  as  bubbles  or  separations  between  soil  clumps  or  rock  fragments.  In  this  condition, 
the  fluid  has  essentially  zero  stiffness  and  will  expand  indefinitely.  Sufficient  volumetric 
compression  can  eventually  collapse  the  void  spaces  such  that  the  pore  fluid  will  again  carry 
compressive  pressures. 
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In  developing  a  model  for  this  behavior,  several  simplifying  assumptions  are  made: 


•  A  single  pore  fluid  model  will  be  used  to  represent  the  pore  water  and  any  gas- 
filled  void  that  may  develop  under  vacuum; 

•  No  attempt  is  made  to  model  the  actual  process  of  gas  void  formation:  the  voids 
are  simply  assumed  to  develop,  without  bounds,  with  expanding  pore  volumes; 

•  The  vapor  pressure  required  to  form  water  vapor  is  neglected; 

•  If  the  pore  fluid  was  partially  saturated  at  the  start  of  the  problem  (see  the  model 
described  in  Section  3.5),  the  fluid  is  assumed  to  return  to  the  initial  partial 
saturation  condition  when  the  voids  created  under  vacuum  have  collapsed. 

In  discussing  this  feature  of  the  water  model,  the  term  "distended"  is  used  to  describe  the 
condition  of  the  pore  fluid  under  vacuum  pressure.  A  "rejoined"  condition  is  said  to  exist  when 
subsequent  compressive  strains  return  the  pore  water  to  the  initial  condition  which  may,  in  fact, 
be  a  partially  saturated  condition. 

Given  the  assumptions,  distention  and  rejoining  of  the  pore  water  can  be  modeled  in  a 
fairly  simple  manner.  A  negative  volume  strain  in  the  pore  fluid,  computed  with  respect  to  the 
initial,  unloaded  condition  indicates  the  onset  of  distention.  While  the  pore  fluid  is  distended, 
the  stiffness  of,  and  pressure  in  the  pore  fluid  is  set  to  zero.  The  fluid  is  considered  to  be 
rejoined  (at  the  initial  conditions)  when  the  pore  water  volume  strain  becomes  compressive 
(positive).  The  total  volume  strain  in  the  pore  fluid  is  computed  from  the  current  porosity  of 
the  material,  but  an  adjustment  is  made  to  account  for  movement  of  water  in  and  out  of  the 
material  as  simulated  in  the  MPDAP  and  MEM  codes. 


3.5  PORE  FLUID  MODEL  FOR  PARTIAL  SATURATION  CONDITIONS 

When  rock  or  soil  is  unsaturated,  compression  of  the  pore  water  and  solid  grains  is 
nearly  insignificant  when  compared  with  the  compression  of  pore  air.  Under  these  conditions, 
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material  behavior  is  governed  mostly  by  the  skeleton  model.  With  sufficient  compression,  the 
pore  air  gets  squeezed  out  and  the  material  becomes  saturated.  Rischbieter,  et.  al.  (1977) 
demonstrated  that  even  a  minute  amount  of  entrapped  air  drastically  alters  the  pore  pressure 
response  in  multiphase  porous  materials.  To  simulate  this  behavior,  the  pore  fluid  model  is 
modified  to  account  for  the  compressibility  of  pore  air  and  converges  to  a  saturated  condition. 
Note  that  this  model  is  invoked  only  when  the  initial  saturation  is  less  than  100%. 


The  compressibility  of  the  air- water  mixture,  C,«,,  is  defined  as: 

C  =  (3.84) 

where  is  the  fluid  pressure.  The  volumetric  strain  in  the  air-water  mixture,  ev..w.  is  the  sum 
of  volume  strain  in  the  air  and  water.  Using  the  definition  of  the  initial  saturation,  it  can  be 
shown  that: 


where 


g  =  (1  -5j  g  +  5  e  (3.85) 

€v,aw  =  volume  strain  of  air-water  mixture 
Cv.a  =  volume  strain  of  air  bubbles 
^v.w  ==  volume  strain  of  water  (from  Equation  3.79) 

So  =  initial  saturation 


From  Equations  3.84  and  3.85  we  can  get  an  expression  for  the  compressibility  of  the  air- water 
mixture: 

c..  =  (1  -SJC,*  S,C,  (3.86) 


Since  the  compressibility  is  the  inverse  of  the  bulk  modulus.  Equation  3.86  can  be  expressed  as: 


1  _  1-^.  , 


(3.87) 


where: 


=  bulk  modulus  of  air-water  mixture 
K,  =  equivalent  bulk  modulus  of  air  bubbles  in  the  fluid 
K„  =  bulk  modulus  of  water  (from  Equation  3.81) 


3-28 


Kim  (1982)  developed  a  formulation  for  the  volume  strain  and  equivalent  bulk  modulus 
of  air  bubbles  in  the  pore  fluid.  This  formulation,  which  was  summarized  by  Kim,  Blouin,  and 
Timian  (1986)  explicitly  modeled  the  dissolution  of  gas  in  liquid.  This  formulation  was 
implemented  and  tested  in  MPDAP.  While  it  is  fundamentally  sound,  its  practical  application 
requires  the  evaluation  of  some  parameters  which  are  nearly  impossible  to  obtain  directly. 
Initially,  the  surface  tension  forces  were  assumed  to  remain  constant.  In  reality,  the  surface 
tension  forces  are  related  to  the  inverse  of  the  bubble  radius.  Thus,  surface  tension  grows 
tremendously  as  the  air  bubbles  are  compressed  and  eventually  collapse.  If  this  effect  is  not 
modeled,  the  effective  modulus  of  the  fluid  is  unrealistically  low.  Since  there  is  no  available 
procedure  for  determination  of  the  bubble  radius,  the  full  theoretical  complexity  of  the 
formulation  is  not  justified. 


As  an  alternative,  a  much  less  complex  model  was  implemented  for  pore  air  based  on  the 
adiabatic  ideal  gas  law; 


It  V  ^ 

ao  ao 


(3.88) 


where; 


TTjo  =  initial  air  pressure  (absolute  pressure) 
T,  =  current  air  pressure  (absolute  pressure) 
=  initial  air  volume 
V,  =  current  air  volume 
7  =  ratio  of  heat  capacity  (Cp/C„) 


For  non-adiabatic  compression  of  air  at  constant  temperature,  7  has  a  value  of  1.0  and  Equation 
(3.88)  reduces  to  Boyle’s  Law. 


The  volume  strain  of  air  can  be  defined  in  terms  of  engineering  strain. 
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(3.89) 


e 


=  1  - 


(K] 

lO 


Substituting  Equation  (3.88)  into  Equation  (3.89),  we  can  express  the  volume  strain  in  the  air 
in  terms  of  air  pressure: 


e 


Vffl 


=  1 


1 


Tl 


a  ) 


Neglecting  the  influence  of  surface  tension. 


(3.90) 


(3.91) 


where: 


T  =  current  pore  water  pressure  (gage  pressure) 
P,  =  reference  atmospheric  pressure 


Substitution  of  Equation  (3.91)  into  Equation  (3.90)  yields 


+  P 


a/ 


The  tangent  bulk  modulus  of  air  bubble  can  be  defined  as 


d  e 


Differentiating  Equation  (3.92)  with  respect  to  r. 


d  e 


_ 


yfi 


VI 
ao 


•i) 


7t  +  P. 


(3.92) 


(3.93) 


(3.94) 


Substitution  of  Equation  (3.94)  into  Equation  (3.93)  yields 
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Y 


lllif  ■  7) 


(3.95) 


3.6  PERMEABILITY  MODEL 


The  finite  element  code  MPDAP  is  capable  of  calculating  the  flow  of  pore  fluid  between 
elements.  The  flow  of  fluid  with  respect  to  the  skeleton  is  controlled  by  Forchheimer’s 
permeability  model  which  is  described  in  detail  in  Section  4.  The  Forchheimer  model  can  be 
expressed  as: 


TT  .  =  -H  W,  +  W.' 


a 


)8 


pA 


(3.96) 


where: 


T  i  =  pore  pressure  gradient 
g  =  acceleration  of  gravity 
Pf  =  mass  density  of  pore  fluid 
p,  =  pore  fluid  viscosity 
a  —  linear  flow  coefficient 
j8  =  quadratic  flow  coefficient 
w  =  apparent  flow  velocity  relative  to  the  skeleton 
U  =  absolute  acceleration  of  pore  fluid 


The  first  term  in  Equation  3.96  is  simply  Darcy’s  law  while  the  velocity  squared  term  was 
apparently  first  proposed  by  Forchheimer  (1901).  The  first  two  terms  represent  the  energy  loss 
due  to  fluid  flow  while  the  last  term  accounts  for  fluid  inertia. 


Equation  3.96  can  also  be  written  in  the  form: 


IT, 


(3.97) 
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where  k'  represents  an  equivalent,  high  pressure  permeability  coefficient  given  by: 


k'  = 


(3.98) 


Hence,  the  flow  of  pore  fluid  in  the  soil  skeleton  is  governed  by  Equations  3.97  and  3.98  and 
the  flow  coefficients  a  and  0  which  can  be  determined  from  laboratory  test  data. 


If  the  soil  skeleton  undergoes  significant  volumetric  strains,  the  permeability  of  the 
skeleton  will  change.  This  is  particularly  true  near  an  explosive  detonation  where  large 
volumetric  compression  restricts  pore  fluid  flow  or,  for  some  materials,  expansion  of  the 
skeleton  yields  greater  fluid  flow.  We  have  implemented  a  simple  empirical  model  to  account 
for  this  behavior  which  varies  the  Darcy  permeability  coefficient  with  changes  in  porosity 


where: 


a  =  a  (3.99) 

o 

Q-o  =  initial  linear  flow  coefficient 
a  =  flow  coefficient  at  a  porosity  of  n 

Cp  =  slope  of  straight  line  relationship  between  porosity  and  the 
logarithm  of  the  permeability 
n  =  current  porosity 
no  =  initial  porosity 


The  form  of  Equation  3.99  is  suggested  by  test  data  reported  by  Blouin  and  Timian  (1986)  and 
shown  in  Figure  3.7.  A  similar  relationship  may  exist  between  porosity  and  the  quadratic 
coefeicient  0  in  the  flow  equation,  but  there  is  insufficient  high  pressure  permeability  data 
available  to  establish  such  a  relationship.  Hence,  0  is  assumed  to  remain  constant  as  the  porous 
skeleton  deforms. 
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Table  3.1.  Parameters  for  the  ARA  Three  Invariant  Material  Model. 

Elastic  Response: 

1/  :  Poisson’s  ratio 
Ky,  ;  elastic  modulus  constant 
n  :  elastic  modulus  exponent 
Ki  :  initial  bulk  modulus  at  low  pressure 
X,  7,  |S  :  unload/reload  bulk  modulus  constants  used  in  ARA  elastic  model 

Compressive  Plastic  Response: 

r  :  elliptical  cap  (compressive  yield  surface)  axis  ratio 
c,  ;  compressive  hardening  constant  for  segment  of  fit 
p,  :  compressive  hardening  exponent  for  segment  of  fit 

Expansive  Plastic  Response: 

T  ;  apparent  tensile  strength  (negative) 

E  :  expansive  yield  constant 
m  :  expansive  yield  exponent 

bi  :  constants  for  segment  "i"  in  a  multilinear  fit  for  defining  Wp 

P,  /  :  constants  in  an  exponential  fit  for  defining  Wp  p^^ 
q  :  expansive  hardening  exponent 

Sj,  r,  :  constants  for  segment  in  a  multilinear  fit  for  defining  tjj 
s,  R,  t  :  constants  in  a  single  line  fit  for  defining  172 
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Table  3.2.  Format  of  input  data  file  specifying  the  ARA  Three  Invariant 
Model  parameters  for  a  given  material. 


Model 

Parameter 

Input  Variable 

Description 

Stress  units 
used  in  fit 

lUNITS 

lUNITS  =  1  Pa 

=  2  dynes/cm- 
=  3  psi 
=  4  psi 
=  5  ksi 
=  6  MPa 
=  7  kg/cm’ 

=  8  Mbar 

T  P 

APEX,  ATMO 

APEX  =  apparent  tensile  strength  (negative) 
ATMO  =  atmospheric  pressure  in  units  of  the 
calculation 

K„,  n,  V 

AKUR,  AN,  APOI 

AKUR  =  elastic  modulus  constant 

AN  =  elastic  modulus  exponent 

APOI  =  elastic  Poisson's  ratio 

r 

AR,  ACRV 

AR  =  elliptical  cap  major/minor  axes  ratio 

ACRV  =  Number  of  pairs  of  c,  and  p, 

describing  compressive  hardening 
function  (max  ACRV  =  4) 

‘•i,  Pi 
(/  =  1  to 
ACRV) 

AACC(0,  AAPC(/) 

AACC(/)  =  compressive  hardening  constant  for 
segment 

AAPC(/)  =  compressive  hardening  exponent 
for  segment  "i" 

> 

AACC(ACRV),  AAPC(ACRV) 

E,  m,  rj2 

AEY,  AMY,  AETAl 

AEY  =  expansive  yield  constant 

AMY  =  expansive  yield  exponent 

AETAl  =  expansive  yield  hardening  constant 

Option  for 
plastic 
potential 
definition 

NETA2 

NETA2  =  1  multi-segment,  linear  fit  of  772  as 
function  of  fp" 

NETA2  =  2  nonlinear  fit  of  772  as  function  of 
fp"  and 

ti 

Si 

IfNETA2  =  1; 

NPTS 

SOCTAC/-),  /•  =  1  to  NPTS 
ATGA(/),  /  =  1  to  NPTS 
ASGA(/),  /  =  1  to  NPTS 

NPTS  =  number  of  pressure  points 

defining  linear  ranges  of  772 
SOCTA(/)  =  pressure  points, 

ATGA(/)  =  772  constant  t^  at  SOCTA(/) 

ASGA(/)  =  772  constant  Sj  at  SOCTA(/) 

I 
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Table  3.2.  Format  of  input  data  file  specifying  the  ARA  Three  Invariant 
Model  parameters  for  a  given  material  (continued). 


Model 

Parameter 

Input  Variable 

Description 

t,  R,  s 

If  NETA2  =  2; 

AATG,  AARG,  AASG 

AATG  =  constant  t  in  nonlinear  tjj  definition 
AARG  =  constant  R  in  nonlinear  tj,  definition 
AASG  =  constant  s  in  nonlinear  tjj  definition 

Option  for 
specifying 

Wp,p*^ 

NWPPK 

NWPPK  =  1  multi-segment,  linear  fit  for 

Wp  pe»k  as  a  function  of  a^^ 

NWPPK  =  2  nonlinear  fit  for  Wp  p^  as  a 
function  of  ao,., 

Si,  bi 

If  NWPPK  =  1 

PDPK(/),  /  =  1  to  4 

WPDPK(/),  /  =  1  to  4 

PDPK(/)  =  pressure  points,  ao<-i 

WPDPK(/)  =  Wp  p,^  at  PDPK(/) 

(values  of  a,  and  b,  implied) 

R  1 

If  NWPPK  =  2 

APBAR,  AL 

APBAR  =  constant  defining  Wp  p^^ 

AL  =  exponent  defining  Wp  p^^ 

1/q 

AALPH 

AALPH  =  inverse  of  expansive  hardening 
exponent  q 

X,  7,  /3,  Kj 

AHLAM,  AHGAM,  AHBET, 
AHBLK 

AHLAM,  AHGAM,  AHBET  =  unload/reload 

constants  in 

ARA  elastic 
model 

AHLAM  =  -1.0  invokes  Lade  and  Nelson 
elastic  model  (AHGAM, 
AHBET  not  used) 

AHBLK  =  initial  bulk  modulus  at  very  low 
pressure 
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Table  3.3  Fluid  compressibility  model  constants 

(see  Section  3.4  for  definitions  of  constants). 


Parameter 

Units 

Fresh  Water 

Sea  Water 

Po 

kg/m^ 

1002.8 

1026 

Cl 

m/s 

1500 

1522 

s, 

- 

2.00 

1.97 

S2 

s/m 

-1.07  X  10-* 

-0.898  X  10-* 

S3 

- 

1.144 

1.123 

Vp. 

m/s 

4000 

4573 

c, 

m/s 

7788 

8653 

'^pi 

MPa 

31,240 

40,600 

I 

I 

I 


Figure  3.1.  Yield  surfaces  used  in  the  ARA  Three  Invariant 
material  model  (after  Dass  and  Merkle,  1986). 
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Figure  3.2.  Fitting  of  ARA  elastic  formulation  to  unload/reload 
hydrostatic  compression  response. 
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I.«i3  (Adiabatic) 
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SECTION  4 

RELATIONSHIP  BETWEEN  MICROSCOPIC 
PROPERTIES  AND  MACROSCOPIC  FLOW  PARAMETERS 


Section  2  of  this  report  describes  analytical  formulations  and  their  computer 
implementations  for  simulating  the  mechanical  behavior  of  a  saturated  porous  material  based  on 
knowledge  of  the  mechanical  properties  of  its  constituent  materials,  i.e.  the  porous  skeleton 
(possibly  soil  or  rock),  the  solid  grains,  and  the  pore  fluid.  An  advanced  three-invariant 
plasticity  model  for  simulating  the  nonlinear  response  of  a  sand  or  rock  skeleton  is  presented  in 
Section  3,  along  with  nonlinear  models  for  solid  grains  and  pore  water  that  are  applicable  to 
very  high  pressures.  Since  the  fluid  pressure  is  determined  by  the  volume  of  pore  space  the 
fluid  occupies,  skeleton  strain  gradients  can  induce  pore  pressure  gradients,  resulting  in  flow  of 
the  pore  fluid.  The  rate  of  pore  fluid  flow  in  response  to  imposed  pressure  gradients  is  a 
function  of  the  propenies  of  both  the  fluid  and  the  porous  (and  permeable)  skeleton.  The  flow 
characteristics  do  not  strongly  influence  the  deformation  response  of  unsaturated  materials. 
However,  in  saturated  porous  materials,  they  can  be  very  important,  particularly  if  the  loading 
is  dynamic. 

This  section  begins  with  a  discussion  of  the  governing  equation  for  fluid  flow  through 
a  porous  material.  With  that  as  a  theoretical  framework,  fluid  flow  test  data  are  presented  for 
a  range  of  porous  materials.  Included  in  the  data  set  are  granular  materials,  both  spherical  and 
angular,  and  porous  limestone.  Finally,  there  is  a  discussion  of  the  relationship  between  the 
permeability  coefficients  and  the  microscopic  geometric  characteristics  of  the  porous  media. 


4.1.  GOVERNING  EQUATIONS  OF  FLOW  THROUGH  POROUS  MEDIA 

Flow  through  porous  materials  has  long  been  studied  by  investigators  in  a  variety  of 
disciplines,  including  hydrology,  petroleum  resci  /oir  engineering,  chemical  engineering  and 
filtering.  In  traditional  civil  engineering,  where  pressure  gradients  and  the  resulting  flow 
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velcKities  are  typically  small,  problems  are  almost  exclusively  formulated  using  the  Darcy  flow 
equation  which  assumes  a  linear  relationship  between  pressure  gradient  and  flow  velocity. 
However,  in  previous  AFOSR-sponsored  work,  Kim,  et.al.  (1988)  reported  that  in  dynamic 
problems  of  interest  such  as  blast-induced  liquefaction,  weapons  effects,  and  possibly 
earthquakes,  pore  pressure  gradients  and  flow  rates  are  of  sufficient  magnitude  that  the  Darcy 
equation  does  not  adequately  describe  the  pore  fluid  response.  To  analytically  treat  the  dynamic 
problems,  Kim,  et.al.  (1988)  adopted  a  more  general  expression  for  pressure  gradient  which 
includes  a  term  that  is  proportional  to  the  square  of  the  fluid  velocity  in  addition  to  the  linear 
(Darcy)  term,  along  with  inertial  terms  related  to  fluid  acceleration.  The  use  of  the  velocity- 
squared  term  was  apparently  first  proposed  by  Forchheimer  (1901),  and  has  been  used  in  various 
forms  by  investigators  over  the  years.  Forchheimer’s  equation,  which  does  not  consider  fluid 
acceleration,  has  the  general  form; 

Tt .  -  aw.  *  bw^  (■^•1) 

where;  rr ;  =  pore  pressure  gradient  in  the  i  direction 

Wj  =  apparent  fluid  velocity  in  the  i  direction 

a,b  =  constants  which  are  functions  of  both  skeleton  and  fluid  properties 

This  formulation,  with  the  addition  of  the  inertial  terms,  was  used  by  Kim,  et.al.  (1988),  and 
again  as  the  framework  for  investigating  the  influence  of  microscopic  skeleton  properties  on  fluid 
flow  by  Blouin,  et.al.  (1990).  While  Forchheimer’s  equation  appears  to  be  virtually  unknown 
in  the  civil  engineering  literature,  its  use  is  widespread  in  chemical  engineering,  particularly  with 
reference  to  porous  metal  filters.  The  constants,  a  and  b,  in  Equation  4-1  are  functions  of  the 
properties  of  both  the  pore  fluid  and  the  porous  medium.  A  more  useful  form  of  the  equation 
will  result  if  it  can  be  expressed  in  terms  of  constants  that  are  dependent  on  either  fluid  or 
porous  medium  properties,  but  not  both.  Green  and  Duwez  (1951)  argue,  based  on  dimensional 
analysis,  that  if  the  factors  controlling  the  pressure  gradient  are  assumed  to  be  fluid  viscosity, 
fluid  density,  a  length  characterizing  the  pore  openings,  and  the  apparent  fluid  velocity,  then 
Forchheimer’s  equation  must  have  the  form: 
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iL  W 

n ,  -  const,  -  ♦  const. 


(4.2) 


& 


where:  /x  =  dynamic  viscosity  of  the  fluid 

Pf  =  mass  density  of  the  fluid 
5  =  unknown  length  related  to  pore  geometry 

If  the  constants  are  combined  with  the  unknown  pore  dimension  term,  5,  and  called  a  and  0, 
the  resulting  expression  has  the  desired  separation  between  quantities  dependent  on  skeleton 
properties  and  quantities  dependent  on  fluid  properties: 


rr  - 


Pf  .2 

+  — i-W- 


(4.3) 


The  relationship  between  the  two  sets  of  coefficients  is: 


(4.4) 


If  the  term  containing  wf  is  neglected  the  result  is  Darcy’s  equation  with  the  constants  in  a 
different  form.  In  some  soil  mechanics  texts,  a  is  called  the  coefficient  of  permeability  or 
absolute  permeability  and  is  often  represented  by  the  symbol,  K.  It  has  units  of  length  squared 
and  the  coefficient  of  the  quadratic  term,  0,  has  units  of  length.  The  relationship  of  a  to  the 
Darcy  permeability,  k,  is  given  by: 


a 


Pf8 


(4.5) 
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where: 


acceleration  of  gravity 


g  = 


In  previous  work,  we  have  referred  to  the  linear  term  as  the  laminar  response  and  the 
quadratic  term  as  turbulent,  an  analogy  to  flow  in  pipes.  However,  Equation  4-3  does  not  have 
the  same  form  as  the  linear  and  turbulent  flow  equations  for  pipes.  In  pipes,  the  quality  of  the 
flow  is  related  to  the  Reynolds  number.  Re,  defined  as: 


Re 


Of 

0 


(4.6) 


where:  V  =  fluid  flow  velocity, 

and  the  other  quantities  are  as  defined  above.  Laminar  flow  occurs  in  pipes  at  low  velocity 
(Reynolds  number),  and  the  pressure  gradient  is  characterized  by  a  linear  dependence  on  flow 
velocity.  At  some  value  of  Reynolds  number,  the  flow  changes  rather  abruptly  from  laminar 
to  turbulent.  If  the  length  quantity  in  the  Reynolds  number  is  taken  as  pipe  diameter,  then  flow 
is  sure  to  be  laminar  if  Reynolds  number  is  less  than  2100,  and  for  practical  purposes,  it  can  be 
assumed  to  be  turbulent  if  Re  is  greater  than  4700.  In  pipes,  the  pressure  gradient  in  the 
turbulent  regime  is  proportional  to  the  square  of  the  fluid  velocity,  and  there  is  no  linear 
dependence  on  velocity.  In  contrast  Forchheimer’s  equation  (4-3),  retains  its  dependence  on 
both  the  first  and  second  powers  of  velocity  over  the  full  range  of  fluid  velocities,  and  the 
relative  importance  of  the  two  terms  varies  according  to  the  details  of  each  specific  problem. 
Green  and  Duwez  (1951)  attributed  the  velocity  squared  term  to  inertial  effects,  related  to  the 
tortuosity  of  the  flow  path.  Bear  (1972)  also  concluded,  based  on  the  work  of  Lindquist  (1933), 
Schneebeli  (1955),  Hubbert  (1956),  and  Scheidegger  (1960),  that  the  appearance  of  the  square 
term  in  Equation  4-3  is  the  effect  of  inertial  forces  associated  with  the  flow  along  the  tortuous 
path  through  the  porous  material.  These  forces  are  proportional  to  fluid  density  and  independent 
of  viscosity,  consistent  with  the  form  of  Equation  4-3.  Thus,  reference  to  turbulent  flow  has 
been  dropped.  Instead,  we  have  adopted  the  nomenclature,  linear  and  quadratic  flow  coefficients 
for  the  quantities  a  and  jS,  respectively. 
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4.2  PERMEABILITY  MEASUREMENTS  ON  GRANXLAR  MATERIALS 


In  order  to  investigate  the  relationship  between  the  macroscopic  fluid  flow  coefficients 
and  the  microscopic  geometric  properties  of  the  porous  medium,  measurements  of  fluid  flow 
were  made  in  a  variety  of  granular  materials.  These  tests  were  conducted  at  a  wide  range  of 
pressures  on  specimens  of  granular  materials  packed  in  a  circular  duct.  :  he  grain  sizes  of  the 
materials  tested  span  a  range  of  more  than  an  order  of  magnitude,  and  both  spherical  and  angular 
grains  shapes  are  represented. 

The  apparatus  used  to  perform  flow  tests  on  the  granular  materials  is  shown  schematically 
in  Figure  4.1.  It  was  fabricated  from  a  standard  hydraulic  cylinder  by  adding  a  sample  tube  and 
the  necessary  additional  ports.  Screens  fastened  to  each  end  of  the  sample  tube  contain  the 
granular  material.  Kerosene  was  selected  as  the  test  fluid  rather  than  water  to  avoid  rusting  the 
.est  apparatus.  The  piston  of  the  hydraulic  cylinder  serves  to  separate  the  test  fluid  from  the 
energy  source  fluid,  which  was  either  hydraulic  fluid  or  high-pressure  nitrogen.  To  expedite  the 
testing,  the  system  was  equipped  with  solenoid  operated  valves  so  that  the  tests  could  be 
conducted  by  a  single  operator  from  the  control  console.  Plumbing  was  also  installed  external 
to  the  test  apparatus  to  allow  the  fluid  to  be  returned  from  the  downstream  side  of  the  test 
specimen  to  the  upstream  reservoir  without  pumping  it  through  the  test  specimen. 

The  tests  were  instrumented  with  two  pressure  transducers  to  measure  fluid  pressure 
upstream  and  downstream  of  the  test  specimen.  In  practice,  the  flow  resistance  between  the 
downstream  end  of  the  specimen  of  granular  material  and  the  low  pressure  reservoir  was  always 
negligible  in  comparison  with  the  flow  resistance  of  the  specimen.  This  was  verified  by  the 
negligible  output  of  the  downstream  pressure  transducer,  and  the  pressure  gradient  was 
calculated  as  the  steady  state  upstream  pressure  divided  by  the  specimen  length.  A  displacement 
transducer  (a  Linear  Variable  Differential  Transformer  or  LVDT)  was  used  to  measure  the 
displacement  of  the  piston.  The  time  derivative  of  the  piston  displacement,  with  the  appropriate 
area  corrections,  is  the  fluid  flow  rate.  The  apparent  fluid  velocity  was  computed  as  the  flow 
rate  divided  by  the  gross  cross  sectional  area  of  the  specimen. 
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Only  the  steady-state  portion  of  each  test  was  considered,  i.e.  where  the  pressure  gradient 
and  flow  rate  were  constant.  Previous  work  experimentally  validated  the  non-steady  state 
portion  of  the  flow  equation  accounting  for  the  acceleration  of  the  pore  fluid  (Kim,  et  al.,  1988). 
Figure  4.2  presents  typical  time  histories  of  upstream  pressure  and  piston  displacement.  The 
data  shown  for  piston  displacement  represent  the  entire  measurement  range  of  the  displacement 
transducer.  The  steady  state  piston  velocity  was  found  by  a  least  squares  fit  to  the  portion  of 
the  LVDT  output  where  the  gage  was  on  scale.  The  steady  state  value  of  upstream  pressure  was 
found  by  averaging  the  pressure  transducer  output  over  that  same  time  window. 

For  each  test  specimen  installed  in  the  permeability  apparatus,  a  series  of  flow  tests  was 
conducted  at  a  range  of  pressure  gradients.  Each  run  yielded  of  values  of  steady  state  pressure 
gradient  and  flow  rate.  Data  from  a  number  of  such  runs  are  plotted  in  Figure  4.3.  When 
Equation  4-1  is  divided  by  w,  the  resulting  expression  plots  as  a  straight  line  on  the  axes  shown 
in  Figure  4.3.  The  slope  of  the  line  is  b,  or  p,//3,  and  its  intercept  is  a.  or  p/a.  For  each 
different  granular  material  tested,  least-squares  fits  to  the  data  as  plotted  in  Figure  4.3  have  been 
used  to  derive  the  linear  and  quadratic  flow  coefficients. 

Figure  4.4  is  a  summary  plot  showing  all  of  the  flow  test  data  acquired  for  granular 
materials  during  the  course  of  this  effort.  The  characteristics  of  the  granular  materials  used  in 
each  test  series  are  presented  along  with  the  resulting  flow  parameters  in  Table  4.1.  The  table 
contains  the  flow  coefficients  both  in  terms  of  a  and  b,  and  in  terms  of  a  and  0.  They  are  most 
conveniently  derived  from  the  test  data  in  terms  of  a  and  b.  However,  for  comparison  with 
other  materials  tested  with  different  fluids,  it  is  desirable  to  express  the  flow  coefficients  as  a 
and  0. 

The  test  data  were  further  analyzed  to  examine  the  influence  of  grain  geometry  on  flow 
characteristics.  Figure  4.5  is  a  log-log  plot  of  the  linear  flow  coefficient,  a,  as  a  function  of 
grain  size  for  materials  consisting  of  spherical  grains.  A  least  squares  fit  to  the  spherical  grain 
data  results  in  the  following  relationship  for  between  grain  diameter  and  the  linear  flow 
coefficient,  a: 
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a  -  8.03  X  10"  ^ 


(4.7) 


where:  a  =  linear  flow  coefficient  with  units  of  in.^ 

d  =  grain  diameter  with  units  of  in. 

The  fit  given  by  Equation  4-7  is  shown  as  a  line  on  Figure  4.5.  The  exponent  on  the  grain 
dimension  in  Equation  4-7  is  nearly  2.  If,  instead  of  an  unconstrained  least  squares  fit,  the  fit 
is  constrained  to  give  the  best  fit  line  with  a  slope  of  2.0,  the  resulting  line  is: 

a  -  0.00104  (‘^•8) 

With  the  exponent  of  d  forced  to  2,  a  and  d^  have  the  same  units  and  the  coefficient,  0.00104, 
is  dimensionless.  Since  a  is  independent  of  fluid  properties,  Equation  4.8  is  applicable  to  flow 
of  any  fluid.  Equation  4-8  is  shown  as  a  dashed  line  in  Figure  4.5.  Over  the  range  of  grain 
sizes  tested,  the  difference  between  this  equation  and  the  unconstrained  least  squares  fit  is 
insignificant. 

Hazen  (1911)  published  the  following  approximate  expression  for  estimating  Darcy’s 
permeability  coefficient,  k: 

k  -  C(£>,o)^  (‘^•9) 


where:  k  =  Darcy  permeability  coefficient  in  cm/s 

C  =  a  constant  in  the  range  1-1.5  cm/(mm^-s) 

Dio  =  size  in  mm  for  10%  pass  on  the  grain  size  curve 

Assuming  for  the  moment  that  C  has  the  value  1  cm/(mm^-s),  then  in  dimensionally 
homogeneous  form,  C  =  10  mm  '-s  *  =  254  in.’‘-s''.  This  can  be  converted  to  an  expression 
for  a  using  Equation  4-5.  Since  the  Hazen  expression  applies  to  water  flow,  the  conversion  is 
made  using  and  the  viscosity  and  density  values  for  water  at  20°C,  fx  =  1.45  *  10'^  Ib-s/in.^  and 
Pf  =  9.35  *  10'^  Ib-s^/in."*,  respectively.  Substituting  Equation  4-9  into  4-5,  and  using  the  value 
of  C  in  inch  units  and  g  =  386.4  in./s^,  results  in  the  following  expression: 
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(4.10) 


a  -  0.001  (D,o)^ 

For  a  material  of  uniform  grain  size,  Djo  is  equal  to  d,  and  this  is  essentially  the  same 
expression  as  Equation  4-8.  The  good  agreement  between  the  data  from  our  high  pressure  flow 
tests  and  the  Hazen  expression  is  an  indication  that  the  current  test  and  analysis  procedures  in 
the  high  pressure  gradient  regime  produce  results  consistent  with  conventional  low  pressure 
engineering  measurement  techniques  for  the  linear  part  of  the  problem.  Returning  to  the  original 
range  of  values  for  C  in  Equation  4-9,  the  constant  in  Equation  4-10,  using  Hazen’s  range  of 
coefficients,  can  have  values  in  the  range  0.001  -  0.0015. 

Figure  4.6  presents  values  of  the  quadratic  flow  coefficient.  /3,  derived  from  flow  tests 
on  spherical  grain  materials.  The  following  expression  represents  a  least  squares  fit  to  the 
quadratic  flow  coefficient  data  as  a  function  of  grain  size: 

p- 0.117  (^•^^) 


where:  0  =  quadratic  flow  coefficient  with  units  of  in. 

d  =  grain  diameter  with  units  of  in. 

If  the  data  are  fit  with  the  exponent  of  d  constrained  to  be  unity,  the  resulting  expression  is: 

p  -  0.067  d  (4.12) 

Similar  to  Equation  4-8,  the  coefficient,  0,  and  the  grain  size,  d,  both  have  units  of  length,  and 
the  constant,  0.067,  is  dimensionless,  making  the  expression  valid  for  0  and  d  in  any  length 
units.  Also,  since  0  is  independent  of  fluid  properties,  it  is  valid  for  any  fluid. 

Green  and  Duwez  (1951)  suggested  that  the  ratio  a/0,  which  has  units  of  length,  is  an 
appropriate  value  to  use  for  5  in  the  computation  of  Reynolds  number  (Equation  4-6).  This 
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would  mean  that  systems  with  the  same  Reynolds  number  would  have  the  same  ratio  of  viscous 
to  inertial  dissipative  forces.  This  definition  of  Reynolds  number  could  be  used  to  design  a  scale 
model  of  a  saturated  porous  system  for  laboratory  experiments.  If  a  scaled  system  could  be 
developed  with  grain  sizes  much  larger  than  the  prototype  materials  of  interest,  then  it  might 
become  practical  to  measure  effects  using  available  instrumentation  that,  in  the  prototype,  would 
be  microscopic.  While  it  may  be  possible  to  relate  the  Reynolds  number  defined  in  this  way  to 
the  onset  of  turbulence  in  the  porous  material,  it  would  not  be  expected  to  occur  at  the  same 
values  of  Reynolds  number  that  apply  to  pipes  since  the  geometry  is  completely  different. 

Tests  were  also  performed  on  materials  with  non-spherical  grains,  including  a  carbonate 
sand  from  Enewetak  Atoll  in  the  Pacific,  and  crushed  marble  and  granite.  Figure  4.7  presents 
a  as  a  function  grain  diameter  for  all  of  the  flow  tests  on  granular  materials,  including  both 
spherical  and  non-spherical  grains.  The  non-spherical  grain  data  scatter  on  both  sides  of  the 
fit  for  spherical  grains.  For  the  linear  flow  coefficient,  a,  the  difference  between  spherical  and 
non-spherical  grains  appears  to  be  within  the  overall  scatter  of  the  data,  indicating  that  its 
dependence  on  grain  shape  is  no  more  significant  that  other  uncertainties  that  remain 
unquantified. 

Figure  4.8  shows  the  relationship  between  quadratic  flow  coefficient,  /?,  and  grain  size 
for  both  spherical  and  non-spherical  grains.  Unlike  the  linear  coefficient,  quadratic  flow 
coefficient  clearly  takes  on  lower  values  for  the  non-spherical  grains  than  in  the  baseline 
spherical  case.  This  means  that  the  resistance  to  flow  through  the  matrix  of  angular  grains 
increases  with  velocity  at  a  greater  rate  than  through  spherical  grains  of  the  same  size.  To  first 
order,  0  for  angular  grains  might  be  taken  as  one  half  the  for  spherical  grains  of  the  same 
size. 

4.3  PERMEABILITY  MEASUREMENTS  ON  SALEM  LIMESTONE 

A  series  of  permeability  tests  was  conducted  on  specimens  of  Salem  (Indiana)  limestone 
under  sponsorship  of  the  Defense  Nuclear  Agency.  Tests  were  conducted  at  a  range  of  pressure 
gradients  to  define  both  the  linear  and  quadratic  flow  coefficients  for  the  material.  To 
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investigate  the  influence  of  rock  skeleton  stress  on  flow  properties,  tests  were  performed  on 
specimens  under  hydrostatic  pressure,  including  pressures  high  enough  to  permanently  crush  the 
rock  cementation. 

For  the  permeability  testing  of  noncohesive  granular  materials  described  in  the  previous 
subsection,  the  test  specimens  were  prepared  by  packing  the  granular  materials  in  a  circular 
duct.  Because  of  the  material’s  granular  nature,  the  packing  process  produced  a  test  specimen 
in  intimate  contact  with  the  wall  of  the  cylinder,  thereby  minimizing  any  error  due  to  flow  along 
the  cylinder  wall.  It  would  be  impractical,  if  not  impossible,  to  achieve  a  close  enough  fit 
between  a  rock  specimen  and  a  steel  cylinder  to  preclude  unwanted  flow  along  the  rock-cylinder 
interface.  Thus,  an  alternative  approach  was  required.  Permeability  testing  of  the  porous  rock 
was  performed  in  a  specially  modified  triaxial  compression  apparatus.  Figure  4.9  is  a  schematic 
depiction  of  the  apparatus  used  to  conduct  permeability  tests  on  porous  limestone.  The  tests 
were  performed  on  cylindrical  rock  specimens  approximately  1-7/8  inch  in  diameter  by  4  inch 
long.  An  impermeable,  flexible  membrane  of  heat  shrinkable  polyolefin  tubing  was  applied  so 
tliat  it  fit  tightly  to  the  rock  and  was  sealed  on  each  end  to  a  perforated  steel  endcap  as  shown 
in  Figure  4.9.  The  specimen,  thus  prepared,  was  installed  in  the  triaxial  compression  apparatus, 
which  was  specially  equipped  with  fluid  flow  passages.  As  illustrated  in  Figure  4.9,  the  top 
and  bottom  endcaps  were  sealed  to  the  piston  and  pressure  vessel  base,  respectively,  with  O- 
rings.  Pressurized  water  was  introduced  to  the  specimen  through  a  port  in  the  base  of  the  vessel 
and  drained  to  the  exterior  of  the  pressure  vessel  through  a  hole  in  the  triaxial  loading  piston. 
In  order  to  assure  that  no  significant  flow  was  allowed  between  the  specimen  and  the  membrane, 
the  triaxial  pressure  vessel  was  pressurized,  and  the  pore  fluid  pressures  were  always  kept  lower 
than  the  confining  pressure  in  the  vessel. 

Flow  generation  and  measurement  techniques  for  the  rock  specimens  were  very  similar 
to  those  used  in  the  tests  on  granular  materials  described  in  Section  4.2.  The  water  was 
pressurized  and  the  flow  rate  measured  with  a  pressure  intensifier  which  is  shown  schematically 
in  Figure  4.9.  The  intensifier  consists  of  two  coupled  hydraulic  cylinders  as  shown  in  the 
figure.  The  cylinders  have  an  area  ratio  of  6.25,  allowing  a  water  pressure  output  of  20,000 
psi  with  an  input  of  hydraulic  fluid  at  3200  psi.  The  intensifier  displacement  was  measured  with 
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an  LVDT  to  determine  the  pore  fluid  flow  rate. 

In  all  cases,  the  specimen  was  first  loaded  hydrostatically  to  a  specified  total  stress  level 
using  a  combination  of  fluid  pressure  and  piston  load.  For  a  simple  hydrostatic  load,  it  would 
not  be  necessary  to  use  the  axial  loading  piston.  However,  for  permeability  testing  axial  loading 
was  applied  with  the  piston  in  order  to  maintain  the  seal  between  the  perforated  top  cap  and  the 
piston  which  was  necessary  to  prevent  flow  of  the  confining  fluid  into  the  pore  space  of  the 
specimen.  Once  the  specified  level  of  hydrostatic  total  stress  on  the  rock  was  reached,  it  was 
held  constant  while  tests  were  performed  at  a  series  of  pore  pressure  gradients.  For  each  of 
those  tests,  the  hydraulic  fluid  flow  to  the  low-pressure  side  of  the  intensifier  was  quickly  raised 
to  the  desired  level  and  held  constant.  When  the  steady  state  condition  (constant  pressure 
gradient  and  flow  rate)  was  achieved,  data  acquisition  was  initiated  to  measure  and  record  the 
pore  fluid  pressure  upstream  of  the  specimen  and  the  fluid  flow  rate.  These  measurements  were 
processed  as  described  in  Section  4.2. 

A  Salem  limestone  specimen  with  porosity  of  0.169  was  the  subject  of  a  series  of  fluid 
flow  tests.  It  was  initially  loaded  to  10  MPa  confining  pressure.  The  first  fluid  flow  test  was 
conducted  with  the  pressure  at  the  upstream  end  of  the  rock  specimen  set  to  approximately 
2.5  MPa.  In  successive  tests,  the  upstream  pore  fluid  pressure  was  incremented  until  it  reached 
a  value  of  about  90%  of  the  confining  pressure.  It  was  not  raised  above  that  level  to  ensure  that 
the  confining  pressure  was  always  higher  than  the  pore  pressure  to  keep  the  membrane  pressed 
firmly  against  the  rock  specimen.  Once  the  pore  fluid  pressure  reached  90%  of  the  current 
confining  pressure,  the  confining  pressure  was  incremented  and  the  test  series  was  repeated  with 
the  pore  pressure  beginning  at  2.5  MPa.  Figure  4.10  presents  data  from  tests  on  a  single 
specimen  at  a  range  of  confining  pressure  and  pore  pressure  gradient  conditions.  In  Figure  4. 10, 
one  line  is  shown  for  each  confining  pressure,  and  each  line  is  identified  by  that  pressure  value. 
The  data  are  plotted,  as  in  Figure  4.3,  such  that  the  intercept  and  slope  of  the  line  are  the  linear 
and  quadratic  flow  coefficients,  a  and  b,  respectively. 

At  pressures  of  70  MPa  and  below,  the  data  points  fall  on  a  straight  line  within  each 
confining  pressure  set,  and  there  is  only  slight  sensitivity  to  changes  in  confining  pressure.  As 
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confining  pressure  is  increased  above  70  MPa,  the  cementation  of  the  rock  skeleton  begins  to 
be  crushed,  thereby  changing  the  geometry  of  the  microscopic  flow  paths  through  the  rock.  This 
results  in  a  more  pronounced  variation  in  flow  coefficients  with  increasing  confinement.  It  is 
also  apparent  from  Figure  4.10  that,  at  the  higher  confining  pressures,  the  lines  are  no  longer 
straight.  Their  slopes  decrease  with  increasing  pore  pressure  gradient.  To  understand  this 
behavior,  it  is  necessary  to  consider  the  stress  conditions  in  the  test  specimen.  The  decrease  in 
permeability  of  the  rock  specimen  with  increasing  confining  pressure  is  caused  by  compaction 
of  the  rock  skeleton  and  the  resulting  reduction  in  the  flow  area  through  its  pore  space.  The 
deformation  of  the  skeleton  is  controlled,  to  first  order  at  least,  by  the  effective  stress,  i.c.  the 
difference  between  confining  pressure  and  pore  pressure.  While  fluid  is  flowing  through  the 
specimen  during  the  test,  there  is  a  pressure  gradient  in  the  pore  fluid  along  the  length  of  the 
specimen.  As  a  result,  there  is  also  a  gradient  in  the  effective  stress  acting  on  the  rock  skeleton, 
with  the  lowest  effective  stress  at  the  upstream  end.  If  the  pore  pressure  at  the  upstream  end 
is  a  small  fraction  of  the  confining  pressure,  then  the  effective  stress  is  approximately  the  same 
as  the  confining  pressure,  and  the  variation  in  effective  stress  along  the  length  of  the  specimen 
is  small  in  comparison  with  the  confining  pressure.  At  the  other  end  of  the  spectrum,  if  the 
upstream  pore  pressure  is  nearly  equal  to  the  confining  pressure,  then  the  effective  stress  at  the 
upstream  end  of  the  specimen  is  nearly  zero,  and  the  variation  in  effective  stress  along  the 
specimen  length  is  very  significant,  being  approximately  equal  to  confining  pressure.  At  low 
confining  pressures,  where  the  rock  skeleton  is  elastic,  and  the  flow  properties  are  relatively 
insensitive  to  changes  in  effective  stress,  the  rock  has  nearly  the  same  flow  properties  along  its 
entire  length  and  the  lines  in  Figure  4.10  appear  nearly  straight.  However,  at  the  higher 
confining  pressures  where  the  skeleton  cementation  is  being  crushed,  i.e.  greater  than  70  MPa, 
the  flow  properties  are  more  sensitive  to  changes  in  effective  stress  and  hence  exhibit  significant 
variation  along  the  specimen  length  in  those  cases  where  the  upstream  pore  pressure  is 
significant  in  comparison  with  the  confining  pressure.  Thus,  for  confining  pressures  above 
approximately  70  MPa,  the  curves  in  Figure  4.10  exhibit  an  initial  linear  range  where  the 
effective  stress  acting  on  the  rock  is  approximately  constant  and  equal  to  the  confining  pressure. 
However,  with  increasing  upstream  pore  pressure,  the  effective  stress  on  the  upstream  portion 
of  the  specimen  is  reduced,  resulting  in  increased  flow  capacity  and  a  downward  curvature  of 
the  curves  in  Figure  4.10. 
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Values  of  the  flow  coefficients  determined  from  the  data  in  Figure  4. 10  are  presented  in 
Table  4.2.  These  values  were  derived  from  the  test  data  using  an  analysis  approach  similar  to 
that  used  for  sands  and  illustrated  by  Figure  4.3.  The  difference  is  that  for  the  limestone,  only 
data  from  tests  in  which  the  upstream  pore  pressure  was  less  that  half  the  confining  pressure 
were  used  in  the  analysis.  As  explained  above,  this  eliminates  those  tests  that  had  large  effective 
stress  gradients  along  the  flow  path  which  tend  to  obscure  the  fundamental  flow  properties  of 
the  rock.  The  coefficients  are  presented  in  the  table  both  in  terms  a  and  b,  as  they  were  derived 
from  the  test  data,  and  in  terms  of  a  and  which  are  properties  only  of  the  porous  rock.  The 
variations  in  the  linear  and  quadratic  flow  coefficients  with  effective  stress  are  illustrated  in 
Figures  4.11  and  4.12,  respectively. 

4.4  RELATIONSHIP  BETWEEN  PERMEABILITY  COEFFICIENTS  AND  SKELETON 

PROPERTIES 

Section  4.2  describes  relationships  between  the  flow  coefficients,  a  and  jS,  and  grain  size 
for  flow  through  porous  materials  consisting  of  spherical  grains  of  approximately  uniform  size, 
and  Figures  4.4  and  4.5  present  the  data  from  which  those  relationships  were  derived.  If  the 
relationships  between  flow  coefficients  and  grain  size  are  derived  based  strictly  on  mathematical 
fits  to  the  finite  sets  of  data  points,  the  results  are  the  empirical  expressions  for  a  and  jS  given 
by  Equations  (4.7)  and  (4.11).  If,  instead,  the  flow  coefficients  are  assumed  proportional  to 
integral  powers  of  grain  size,  as  suggested  by  dimensional  analysis,  the  resulting  relationships 
are  those  given  by  Equations  (4-8)  and  (4-12).  In  the  latter  set  of  equations,  the  numerical 
constants  are  dimensionless. 

The  grain  size  was  a  convenient  independent  variable  for  analysis  of  granular  materials. 
However,  as  attention  is  turned  coward  identifying  a  more  fundamental  relationships  between 
microscopic  geometric  properties  and  permeability  for  various  types  of  materials,  it  is  no  longer 
suitable.  A  conceptual  model  that  is  widely  used  for  the  study  of  flow  through  porous  media 
is  a  matrix  of  openings,  or  pores,  interconnected  by  smaller  channels,  or  throats.  For  spherical 
grains,  the  throats  have  dimensions  of  the  order  of  the  grain  size.  The  effective  diameter  of  a 
throat  was  estimated  as  the  diameter  of  the  circle  whose  area  is  equivalent  to  the  actual  area  of 
the  throat.  If  the  spherical  grains  are  assumed  to  be  close  packed,  the  smallest  throats  will  be 
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formed  by  the  three  tightly  packed  spheres,  and  the  equivalent  diameter  of  the  throat  is  0.23  d. 
Assuming  that  the  throat  diameter,  t,  is  equal  to  0.23  d,  equations  4.8  and  4. 12  can  be  rewritten 
in  terms  of  t  as: 

a  -  0.0197 

p  -  0.509  r  (4.14) 

It  is  of  interest  to  see  how  well  Equations  (4-13)  and  (4-14),  which  were  derived  from 
test  data  from  granular  materials,  predict  the  permeability  properties  of  the  porous  limestone. 
Wardlaw,  et.al  (1987)  reported  the  throat  dimensions  of  two  different  rock  types,  including 
Indiana  limestone,  derived  from  Wood’s  metal  porosimetry.  In  this  technique,  molten  Wood’s 
metal,  a  low  melting  point  (70  C)  alloy  of  bismuth,  is  used  to  saturate  the  pore  space  of  the 
rock.  The  Wood’s  metal  saturated  specimen  is  then  sectioned  and  polished,  and  the  rock 
skeleton  is  dissolved  away  with  acid  leaving  the  metal  matrix  in  the  shape  of  the  pore  space  of 
the  rock.  The  degree  of  saturation  achieved  is  a  function  of  the  pressure  applied  to  force  the 
molten  metal  into  the  rock.  The  highest  degree  of  saturation  reported  by  Wardlaw,  et.al.  was 
94%,  a  result  of  the  application  of  1700  psi  saturation  pressure.  Measurements  were  made  of 
the  resulting  metal  matrix  on  enlarged  scanning  electron  micrographs. 

In  480  observations,  Wardlaw,  et.al.  reported  throats  ranging  from  2  x  to  8  x  lO"* 
inch  with  an  average  of  1  x  10'’*.  If  it  is  assumed  that  the  flow  is  controlled  by  the  smallest 
throat  size,  and  a  and  (3  are  computed  using  a  throat  size,  t,  of  2  x  10'^  inch  in  Equations  4-13 
and  4-14,  the  resulting  permeability  coefficients  are  a  =  7.88  x  10  '^  in^  and  /3  =  1.02  x  10'^ 
in.  The  linear  permeability  coefficient  value  of  7.88  x  10  '^  in^  is  in  remarkably  good  agreement 
with  the  measured  value  of  1.26  x  10  "  in^.  Admittedly,  the  choice  of  the  minimum  measured 
pore  size  as  input  to  Equation  4-13  was  somewhat  arbitrary.  However,  we  believe  that  this 
result  suggests  a  strong  dependence  of  a  on  the  throat  dimension  as  would  be  expected  for  a 
viscosity  related  term.  In  contrast,  the  value  of  the  computed  quadratic  flow  coefficient,  /3  = 
1.02  X  10'^  inch,  differs  by  about  four  orders  of  magnitude  from  the  measured  value  (at  low 
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effective  stress)  of  1.2  x  lO'*  inch.  This  limited  evidence  appears  to  suggest  that  /3  has  a 
stronger  dependence  on  some  characteristic  of  the  porous  skeleton  other  than  throat  diameter. 

This  finding  is  consistent  with  the  hypothesis  that  the  quadratic  flow  coefficient  is 
dependent  on  the  tortuosity  of  the  flow  path  through  the  porous  material.  For  the  granular 
materials,  an  approximate  linear  correlation  was  found  between  the  grain  size  and  the  quadratic 
flow  coefficient,  iS.  However,  if  tortuosity  can  be  characterized  by  a  quantity  with  units  of 
length,  such  as  curvature  (1/length),  then  it  scales  with  grain  size.  Thus,  if  is  proportional 
to  inverse  curvature,  then  for  granular  materials,  it  will  also  be  proportional  to  grain  size,  even 
if  it  is  not  the  size  of  the  throats  that  determine  /3.  Since  the  pore  geometry  of  the  rock  has 
completely  different  shape  characteristics  than  the  granular  materials,  it  would  not  be  expected 
to  scale  with  throat  size  in  the  same  way  as  the  sands,  and  hence  the  expression  relating  13  to 
throat  dimension  for  granular  materials  would  not  be  applicable  to  the  rock. 

4.5  RELATIVE  IMPORTANCE  OF  LINEAR  ANT)  QUADRATIC  TERMS 

The  existence  of  the  quadratic  term  in  the  Forchheimer  expression  for  relating  pressure 
gradient  to  fluid  flow  velocity  can  readily  be  demonstrated  in  the  laboratory.  However,  it  is 
apparently  not  widely  known  or  used  in  civil  engineering.  This  leads  to  the  question  of  where 
it  is  significant.  Figure  4.13  represents  an  attempt  to  shed  some  light  on  that  question.  The 
figure  shows  the  relative  contributions  of  the  linear  and  quadratic  terms  to  the  pressure  gradient 
for  two  sands  and  a  porous  rock  under  a  range  of  flow  conditions.  Table  4.3  presents  the 
properties  of  the  materials  used  to  compute  the  data  points  in  Figure  4.13.  The  two  sands  were 
assumed  to  have  grains  that  would  just  pass  U.S.  Standard  No.  10  and  No.  100  sieves.  Flow 
coefficients  for  the  sands  were  computed  using  Equations  4-13  and  4-14.  The  limestone  flow 
coefficients  were  taken  from  the  low  flow  end  of  the  test  data  presented  in  Table  4.2.  In  all 
cases,  the  flow  coefficients  were  assumed  to  remain  constant  over  the  entire  range  of  pressure 
gradients  considered.  The  fluid  was  assumed  to  be  water  at  20°  C,  and  its  properties  are 
included  in  Table  4.3.  The  variable  on  the  vertical  axis  in  Figure  4.13  is  the  ratio  of  the 
quadratic  to  linear  terms  in  Equation  4-1.  The  figure  illustrates  the  range  of  pressure  gradients 
at  which  the  quadratic  term  becomes  important  for  the  various  materials. 
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A  horizontal  line  is  shown  on  Figure  4.13  at  a  ratio  of  0.10.  At  this  point,  the 
contribution  of  the  quadratic  term  is  10%  of  the  linear  term  contribution.  Taking  that  as  a 
criterion  for  significance,  flow  conditions  above  the  line  on  the  figure  are  significantly  influenced 
by  the  quadratic  term  which  is  missing  from  the  traditional  Darcy  permeability  formulation. 
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Table  4.1.  Summary  of  flo^v  test  data  for  samples  of  uniform  granular  materials  tested  with  kerosene. 


lote:  for  kerosene  at 


Table  4.2.  Permeability  Test  Data  for  Salem  Limestone 


Effective 

Stress 

(psi) 

a 

(lb-s/in^4) 

b 

(lb-s^2in^5) 

Flow  Coefficients 
alnha  beta 

{in^2)  (in) 

1450 

1.26E+C4 

n 

8.79E-H04 

1.15E-11 

1.06E-09 

2175 

1.33E+04 

8.32E-I-04 

1.09E-11 

1.12E-09 

2900 

1.47E-t-04 

7.37E-i-04 

IHlf 

1.27E-09 

3625 

1.54E+04 

■'.57E+04 

9.41E-12 

1.27E-09 

4350 

1.71E+04 

7.13E-I-04 

8.49E-12 

1.31E-09 

5800 

1.86E+04 

1.02E +05 

9.15E-10 

7250 

2.26E+04 

1.21E+05 

7.72E-10 

10150 

3.09E-f-04 

2.47E+05 

4.71E-12 

3.78E-10 

13050 

4.95E+04 

6.91E+05 

2.94E-12 

1.35E+-10 

15950 

9.13E+04 

1.49E+06 

1.59E-12 

6.30E-11 

18850 

1.41E+05 

5.15E+06 

1.03E-12 

i.8:e-ii 

21750 

2.88E+05 

l.llE+07 

5.04E-13 

8.39E-12 

Note:  for  water  at  20°c, 


/i  =  1.45  X  10'^  Ib’s/in^ 


p  =  9.35  X  10-^  Ib-sViiT* 
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Table  4,3.  Properties  of  materials  used  to  create  Figure  4.13. 


a 

(in’) 

0 

(in.) 

a 

lb  •s/in'* 

Sand  No.  10 

6.20  X  10" 

0.0053 

0.0235 

0.018 

Sand  No.  100 

3.48  X  10'" 

0.00040 

4.17 

0.236 

Salem  Limestone 

1.26  X  10  “ 

1.2  X  10-’ 

12161 

71250 

Note:  Water  at  20°c 

=  1.45  X  10'^  lb  •  s/in- 
Pf=  9.35  X  10"*  Ib-inVs^ 
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Figure  4.1.  Schematic  section  view  of  apparatus  used  to 
measure  flow  through  granular  materials. 
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Quadraiic  Flow  Coefficient 


Linear  Flow  Coefficient, 


Quadratic  Flow  Coefficient 


Figure  4.9.  Schematic  of  flow  test  apparatus  for  rock  specimens. 
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1.  Variation  in  linear  flow  coefficient,  ot 
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.12.  Variation  in  quadratic  flow  coefficient,  6, 
with  effective  stress  in  Salem  limestone. 
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Figure  4.13. 


Relative  importance  of  the  quadratic  and  linear 
flow  terms  in  the  Forchheimer  equation  for  two 
sands  and  a  porous  rock. 


SECTION  5 


ANALYSIS  OF  QUASI-STATIC  UNDRAINED,  TRIAXIAL 
COMPRESSION  RESPONSE  OF  SALEM  LIMESTONE 

Applied  Research  Associates  conducted  a  series  of  laboratory  material  property  tests  on 
drained  and  undrained  samples  of  porous  limestone  from  the  Salem  formation  in  Indiana.  This 
work  was  performed  under  sponsorship  of  the  Defense  Nuclear  Agency  (DNA)  with  the  goal 
of  providing  data  to  formulate  and  validate  two-phase  models  of  the  limestone  response.  The 
results  of  those  tests,  which  are  reported  by  Chitty  &  Blouin  (1993)  and  summarized  by  Kim, 
Blouin,  Chitty  and  Merkle  (1988),  provide  an  excellent  opportunity  to  exercise  the  analytical 
formulations  developed  under  AFOSR  sponsorship. 

In  the  DNA  effort,  drained  tests  were  performed  to  define  the  rock  skeleton  properties. 
Based  on  those  skeleton  properties,  two-phase  analytical  models  were  derived  to  characterize  the 
saturated  undrained  behavior  of  the  material.  In  order  to  validate  the  two-phase  analytical 
models,  tests  were  also  performed  in  the  saturated  undrained  condition,  including  undrained 
triaxial  compression  tests.  The  undrained  triaxial  compression  test  results  contained  apparent 
inconsistencies  that  motivated  further  study.  Specifically,  the  test  data  showed  that  pore  pressure 
remained  essentially  constant  while  the  sample  underwent  several  percent  of  volumetric 
expansion.  In  an  effort  to  better  understand  the  pore  pressure  and  volumetric  strain  behavior 
of  these  tests,  ARA  performed  a  supplemental  study,  sponsored  by  the  U.S.  Army  Engineers 
Waterways  Experiment  Station  (WES),  in  which  additional  instruments  were  used  to  better 
define  the  deformed  shape  of  the  test  specimens  (Chitty  &  Blouin,  1990). 

The  objectives  of  the  analysis  described  here  were  to  exercise  and  validate  the  material 
modeling  capabilities  of  MEM  as  developed  for  this  project,  and  to  gain  a  fundamental 
understanding  of  the  undrained  triaxial  compression  test  and  its  applicability  to  undrained 
material  property  characterization. 

This  section  presents  an  analysis  of  the  observed  response  of  Salem  limestone  in  a 
conventional,  undrained  triaxial  compression  test.  First,  the  appropriate  model  parameters  used 
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in  numeiical  simulations  of  the  triaxial  test  are  presented.  Then,  the  assumptions  and  data 
analysis  procedures  that  are  conventionally  applied  to  triaxial  compression  tests  of  rock  are 
described  together  with  the  results  of  undrained  tests  on  Salem  limestone  including  special 
deformation  measurements  that  were  acquired.  Finally,  the  results  of  a  finite  element  simulation 
a  the  laboratory  test  are  described  followed  by  a  discussion  of  the  conclusions  resulting  from  this 
investigation.  The  results  presented  in  this  section  were  summarized  in  a  paper  by  Chitty,  et. 
al,  (1991). 

5.1  MODEL  PARAMETER  DETERMINATION  FOR  SALEM  LIMESTONE 

Salem  limestone  is  a  porous  rock  of  exceptional  uniformity.  The  samples  used  in  the 
undrained  triaxial  tests  had  an  average  porosity  of  0.128.  Table  5.1  summarizes  the  physical 
properties  of  the  material  used  in  the  numerical  simulations. 

The  solid  grains  of  the  limestone  were  modeled  using  the  parameters  given  in  Table  5.2 
together  with  the  numerical  model  described  in  Section  3.3.  The  limestone  was  saturated  with 
fresh  water  in  the  simulations;  the  partial  saturation  model  described  in  Section  3.5  was  not 
used.  The  pore  water  compressibility  was  modeled  as  nonlinear  corresponding  to  the  model 
described  in  Section  3.4  and  the  parameters  given  in  Table  3.3.  The  initial  bulk  modulus  of  the 
pore  water  is  2220  MPa.  Given  that  the  triaxial  tests  were  quasi-static  with  very  low  pore  water 
velocities  expected,  the  Forchheimer  permeability  model  was  not  necessary.  A  representative, 
constant  Darcy  permeability  of  4.70  x  10"’  m/s  was  used  in  all  of  the  numerical  simulations. 
This  value  is  higher  by  approxiately  a  factor  of  6  than  the  test  data  presented  in  Section  4.  The 
value  used  in  the  analysis  was  based  on  a  very  preliminary  laboratory  test  which  was  shown  to 
be  significantly  in  error  as  a  result  of  the  subsequent,  more  precise  tests  described  in  Section  4. 
This  error  does  not  significantly  affect  the  conclusions  of  this  anlaysis  because  the  permeabilities 
in  both  the  triaxial  compression  test  and  the  MEM  calculation  were  high  enough  that  any  pore 
pressure  gradients  equilibrated  essentially  instantly  on  the  time  scale  of  the  test. 
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5.1.1  Fitting  of  ARA  Three  Invariant  Skeleton  Model 

The  determination  of  the  parameters  for  the  ARA  3-1  model  to  represent  the  skeleton 
response  of  Salem  limestone  is  described  in  this  subsection.  The  model  was  fit  to  the  data  from 
the  drained  hydrostatic  and  triaxial  compression  tests  shown  in  Figures  5. 1  and  5.2.  The  model 
parameters  are  summarized  in  Table  5.3. 

Using  the  measured  response  of  Salem  limestone  in  the  elastic  regime  during  an 
unconfmed  compression  test,  a  Poisson’s  ratio  of  0.25  was  determined.  The  initial  bulk  modulus 
measured  in  the  linear  ponion  of  the  drained  hydrostatic  loading  in  Figure  5.1  was  24,000  MPa. 
This  corresponds  to  an  initial  Young’s  modulus  of  36,000  MPa  for  a  Poisson’s  ratio  of  0.25. 
Also,  from  the  cycled  load/unload  hydrostatic  compression  data  in  Figure  5.1,  values  of  the 
initial  unload  bulk  moduli  at  various  maximum  pressures  were  determined,  converted  to  Young’s 
moduli  according  to  Equation  3.11,  and  plotted  as  shown  in  Figure  5.3.  A  least  squares  linear 
regression  was  then  used  to  fit  the  line  in  log-log  space  as  shown  in  Figure  5.3.  This  is 
equivalent  to  the  procedure  outlined  in  Section  3.2.3  and  gives  the  elastic  parameters  n  = 
0.2286  and  =  71,560  as  given  in  Table  5.3. 

The  option  for  the  ARA  elastic  model  was  selected  to  match  the  unloading  response  of 
the  limestone  as  closely  as  possible.  The  unloading  cycle  from  200  MPa  during  the  hydrostatic 
compression  test  (Figure  5.1)  was  used  to  obtain  the  parameters: 

X  =  0.644 
7  =  0.893 
/3  =  0.487 

These  parameters  are  defined  in  Equations  3.4,  3.6  and  3.7.  This  completes  the  definition  of 
the  elastic  response  for  the  skeleton  model. 

The  hydrostatic  compression  data  from  Figure  5. 1  are  also  used  to  define  the  compressive 
plastic  response  as  outlined  in  Section  3.2.4.  The  ratio  of  the  major  to  minor  axes  of  the 
elliptical  compressive  yield  surface  (r  in  Equation  3.15  and  Figure  3.1)  was  selected  to  be  r  = 
2.0  based  on  previous  experience.  The  parameter  fitting  routines  in  MEM  were  then  used  to 
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compute  the  plastic  compressive  work  from  the  plastic  volume  strains  according  to  Equation  3.17 
in  the  manner  described  in  Section  3.2.4.  The  normalized  compressive  plastic  work  was  plotted 
as  a  function  of  pressure  as  shown  in  Figure  5.4.  Four  line  segments  were  then  fit  to  this  data 
defining  four  sets  of  the  parameters  C;  and  p;  as  given  in  Equation  3.20.  The  resulting 
parameters  are  given  in  Table  5.3  and  complete  the  definition  of  the  compressive  plastic 
response. 

The  first  step  in  determining  the  expansive  plastic  response  parameters  is  to  define  the 
ultimate  failure  surface  and  the  corresponding  peak  expansive  plastic  work  as  functions  of  mean 
stress.  The  failure  envelope  is  fit  to  the  ultimate  shear  strengths  from  the  triaxial  compression 
test  data  of  Figure  5.2a  in  the  conventional  manner.  Below  the  brittle  ductile  transition,  the 
ultimate  octahedral  shear  stress  is  used.  In  the  ductile  regime,  the  octaliedral  shear  stress  at 
15%  axial  strain  is  used. 

The  expansive  plastic  response  is  fit  as  outlined  in  Section  3.2.5.  The  failure  envelope 
is  fit  by  plotting  the  failure  points  as  shown  in  Figure  5.5  (note  that  Uec,  =  ^oci  •  T  Figure 
5.5).  Then,  according  to  Equation  3.40,  m  is  determined  to  be  4.203  x  10‘®  from  the  ratio  of 
the  slope  to  the  intercept  of  the  straight  line  fit  to  the  data  in  Figure  5.5.  Then,  from  available 
test  data  on  Salem  limestone,  a  value  of  =  1.25  was  determined  which  gives  E  =  0.1111 
according  to  Equation  3.23.  Then,  substitution  of  m  and  E  into  the  intercept  value  from  Figure 
5.5  gives  ?7j  =  0.4197  (refer  to  Equation  3.40).  The  tensile  strength  of  the  skeleton  was 
measured  to  be  T  =  32  MPa.  These  parameters  are  summarized  in  Table  5.3. 

The  next  step  in  the  fitting  process  is  to  define  the  function  for  the  peak  plastic  expansive 
work  using  a  multi-segment  fit  (Equation  3.28)  or  using  an  exponential  function  (Equation  3.29). 
We  elected  to  use  the  multi-segment  function  to  permit  as  accurate  a  fit  as  possible.  The  peak 
plastic  work  was  computed  using  the  MEM  fitting  routines  to  integrate  the  plastic  expansive 
strains  up  to  the  failure  surface  according  to  Equation  3.27.  The  peak  plastic  work  from  each 
test  was  then  plotted  as  shown  in  Figure  5.6  and  three  straight  line  segments  were  fit  using  least- 
squares  linear  regression.  The  fit  is  specified  in  terms  of  the  four  pressure-peak  plastic  work 
data  pairs  given  in  Table  5.3  which  correspond  to  the  following  parameters  for  Equation  3.28: 
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Segment 

1 

0.0037 

0.118 

2 

0.10 

-8.63 

3 

0.098 

-8.33 

The  next  step  is  to  define  the  expansive  yield  hardening  function  given  by  Equation  3.25. 
Since  the  coefficient  171  and  the  peak  plastic  expansive  work  have  been  previously  defined,  only 
the  exponent  q  is  needed  to  define  the  hardening  function.  In  accordance  with  Equation  3.42, 
the  triaxial  test  data  is  plowed  as  shown  in  Figure  5.7  using  the  MEM  fitting  routines.  Then, 
q  is  computed  at  a  value  of  fp'VTj,  =  0.80  and  plotted  as  a  function  of  mean  stress  in  Figure  5.8. 
A  representative  value  for  q  is  the  average  value  of  about  0.4.  This  yields  the  model 
relationship  plotted  in  Figure  5.9  which  is  a  good  approximation  of  the  test  data  plotted  in 
Figure  5.7. 

The  final  step  in  fitting  the  ARA  3-1  model  is  defining  the  expansive  flow  rule  and 
expansive  plastic  potential  function  as  defined  in  Equations  3.34  and  3.35.  Of  the  parameters 
in  these  functions,  only  the  coefficient  772  is  yet  to  be  determined.  For  Salem  limestone,  we 
elected  to  define  173  as  a  multi-segment  function  of  the  expansive  yield  hardening  (Equation 
3.36).  As  outlined  in  Section  3.2.5,  772  is  computed  from  the  expansive  plastic  strains  using 
Equations  3.44  through  3.46  for  each  triaxial  test. 

Expansive  plastic  volumetric  strains  are  plotted  in  Figures  5.10a  through  5.10d  as  a 
function  of  shear  strain  for  drained  triaxial  test  data  at  confining  stresses  of  50,  100,  200,  and 
400  MPa.  These  plots  are  generated  using  the  fitting  routines  in  MEM.  In  the  case  of  low 
confining  stresses  (50  and  100  MPa),  the  expansive  plastic  volumetric  strains  are  initially 
compressive  and  then  continuously  dilative  as  the  shear  strain  increases.  In  the  case  of  high 
confining  stresses  (200  and  400  MPa),  however,  the  expansive  plastic  volumetric  strains  show 
continuous  dilation  with  no  initial  compression.  To  obtain  representative  average  slopes  for  the 
computation  of  772,  the  expansive  plastic  strain  data  have  been  fit  with  segments  described  by 
quadratic  least  squares  fits  as  shown  by  the  solid  lines  in  Figure  5.10.  The  plastic  expansive 
strains  are  also  compared  to  the  total  volume  strain  in  Figure  5.10. 
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The  plastic  potential  parameter  tjj  (Equation  3.44)  is  then  computed  and  plotted  as  a 
function  of  fp"  (Equation  3.25)  in  Figure  5.11.  The  expansive  plastic  potential  coefficients  tj  and 
Sj  (Equation  3.36)  are  then  found  by  a  least  squares  fit  to  the  data  as  shown  in  Figure  5.11.  The 
computed  values  of  t^  and  Sj  are  plotted  as  a  function  of  the  average  mean  stresses  in  Figures 
5. 12  and  5.13,  respectively.  Both  t;  and  S;  change  rapidly  at  mean  stresses  near  200  MPa  where 
the  limestone  undergoes  a  change  of  failure  mode  from  brittle  to  ductile.  Fits  to  these  data  are 
approximated  as  straight  line  segments  with  the  plastic  potential  coefficients  given  in  Table  5.3. 
This  completes  the  specification  of  all  the  parameters  required  for  representing  the  limestone 
skeleton  with  the  ARA  3-1  model. 

5.1.2  Prediction  of  Drained  Limestone  Response 

In  this  subsection  the  fit  of  the  ARA  3-1  model  is  verified  by  comparisons  between 
numerical  simulations  and  laboratory  tests.  First,  MEM  is  used  to  compute  a  set  of  synthetic 
drained  hydrostatic  and  triaxial  compression  test  data  for  comparison  to  the  data  from  which  the 
material  parameters  are  extracted.  Next,  the  drained  uniaxial  strain  response  of  Salem  limestone 
is  predicted  and  compared  to  drained  uniaxial  strain  test  data.  This  last  comparison  serves  as 
a  check  of  the  model  and  material  parameters  against  a  set  of  test  data  obtained  along  a  strain 
path  not  used  in  developing  the  model  parameters. 

The  drained  hydrostatic  compression  response  computed  by  MEM  is  shown  in 
Figure  5.14  along  with  the  laboratory  data.  The  calculation  is  in  excellent  agreement  with  the 
laboratory  data.  It  should  be  noted  that  for  this  calculation,  the  expansive  part  of  yield  surface 
in  the  model  is  not  active  since  material  is  in  pure  isotropic  compression. 

A  series  of  four  MEM  calculations  were  performed  to  predict  the  drained  triaxial 
compression  response  at  confining  stresses  of  50,  100,  200  and  400  MPa.  The  computed  results 
are  presented  in  Figure  5.15  through  5.18  along  with  the  corresponding  laboratory  test  data  for 
the  corresponding  confining  stresses.  For  each  confining  stress,  both  stress  differences  and 
volumetric  strains  are  plotted  as  functions  of  axial  strain.  Overall,  the  MEM  calculations  agree 
well  with  the  laboratory  drained  triaxial  compression  test  results.  In  the  case  of  200  MPa 
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confining  stress,  the  MEM  calculations  shown  in  Figure  5.17b  overpredict  the  peak  volumetric 
strain  by  about  20% .  This  overprediction  is  probably  due  to  the  coarse  linear  approximation  of 
the  expansive  potential  parameter  7/2  at  the  transition  pressures  between  150  MPa  and  280  MPa. 

The  last  MEM  computation  in  this  section  is  the  drained  uniaxial  strain  response 
prediction  shown  in  Figure  51.9  and  compared  with  test  data.  The  computed  axial  stresses  are 
slightly  underpredicted  in  the  elastic  region  and  are  slightly  overpredicted  beyond  the  elastic 
limit.  Considering  the  ARA  3-1  material  parameters  are  solely  determined  fr~  ..  the  drained 
hydrostatic  and  triaxial  compression  tests,  such  close  agreement  indicates  strong  potential 
applicability  of  the  model  to  other  more  complex  loading  problems. 


5.2  CONVENTIONAL  TRIAXIAL  COMPRESSION  TESTING 

In  a  triaxial  compression  test,  a  right  circular  cylinder  of  material  is  first  loaded 
isotropically  to  a  predetermined  value  of  confining  pressure.  That  confining  pressure  is  then 
held  constant  while  a  compressive  strain  is  imposed  along  the  axis  of  the  cylinder.  The  axial 
strain  is  applied  to  the  sample  by  a  hardened  steel  piston  and  the  sample  is  supported  at  the 
opposite  end  by  a  similar  hardened  steel  platen.  Several  triaxial  compression  tests  at  different 
confining  pressures  are  typically  used  to  define  a  strength  envelope  for  a  material. 

Since  the  elastic  modulus  of  the  steel  end  caps  is  greater  than  that  of  a  typical 
geotechnical  sample,  there  is  a  mismatch  in  radial  expansion  between  the  sample  and  the  mating 
steel  surfaces,  with  the  sample  tending  to  expand  more  than  the  steel.  If  the  interface  were 
perfectly  frictionless,  this  mismatch  in  radial  deformation  would  not  be  significant  since  one 
surface  would  simply  slide  relative  to  the  other  leaving  a  state  of  pure  normal  stress  acting 
across  the  interface.  However,  because  there  is  friction,  the  unequal  deformation  imposes  shear 
stresses  on  the  ends  of  the  sample  in  addition  to  the  desired  axial  stress.  These  shear  stresses 
are  directed  approximately  toward  the  center  of  the  sample  and  have  the  effect  of  increasing  the 
effective  confinement  at  the  ends  of  the  sample.  As  a  direct  result  of  this  apparent  elevated 
confinement  near  the  sample  ends,  a  typical  triaxial  compression  sample  has  a  nonuniform 
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distribution  of  radial  strain  along  its  axis,  with  substantially  more  radial  expansion  near  the  mid¬ 
height  of  the  sample  than  at  its  ends,  i.e.  it  becomes  barrel  shaped.  This  barreling  effect  can 
be  sharply  reduced  in  soils  testing  by  lubricating  the  surfaces  that  mate  with  the  sample. 
However,  in  rock  testing,  where  tests  are  typically  conducted  at  significantly  higher  pressures, 
lubrication  was  not  in  general  use. 

In  order  to  account  for  nonuniform  strains  in  test  data  interpretation,  it  is  typically 
assumed  that  the  test  section  of  interest  consists  of  the  central  third  of  the  test  cylinder.  In  that 
region  the  radial  stress  is  equal  to  the  confining  pressure  and  the  influence  of  end  constraint  is 
considered  negligible.  Under  this  assumption,  the  radial  strain  measured  at  sample  mid-height, 
where  radial  expansion  is  greatest,  is  considered  representative  of  the  deformation  under  the 
nominal  stress  conditions.  The  strains  in  the  radial  and  tangential  directions  are  assumed  to  be 
equal  as  they  would  be  if  the  loading  were  applied  by  a  perfectly  uniform  normal  traction  on  the 
ends  and  along  the  sides.  Further,  since  axial  strain  is  typically  derived  from  a  deformation 
measurement  over  the  entire  length  of  the  sample,  the  axial  strain  is  assumed  to  be  uniform 
throughout  the  sample. 

While  there  may  be  justification  for  these  assumptions  in  dry  or  drained  testing,  they 
immediately  lead  to  problems  when  applied  to  an  undrained  test.  The  strength  measured  in  an 
undrained  triaxial  compression  test  is  a  function  of  the  mean  effective  stress,  which  is  equal  to 
the  difference  between  total  mean  stress  and  the  pore  pressure.  The  pore  pressure  depends  on 
the  change  in  volume  of  the  pore  spaces,  which  is  closely  related  to  volume  strain  Thus,  the 
nonuniform  radial  strains  and  the  corresponding  nonuniform  volume  strains  result  in  nonuniform 
development  of  pore  pressures  in  the  sample.  If  the  sample  is  highly  impermeable,  a  pore 
pressure  gradient  will  be  induced.  If,  on  the  other  hand,  the  sample  is  sufficiently  permeable 
to  allow  the  pore  pressure  to  equilibrate  as  the  test  progresses,  then  the  pore  pressure  depends 
on  the  volume  strain  averaged  over  the  entire  sample  and  not  at  any  one  position.  In  either  case, 
as  illustrated  by  Figure  5.20,  the  volume  strain  measured  at  the  center  of  the  sample  does  not 
correspond  to  the  pore  pressure  measured  at  one  end. 
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5.3  SPECIAL  TEST  TO  STUDY  PORE  PRESSURE  AND  VOLUME  STRAIN 


In  an  effort  to  understand  the  pore  pressures  measured  in  undrained  triaxial  compression 
tests,  ARA  conducted  a  small  WES-sponsored  test  program  to  measure  additional  details  of 
sample  deformation.  Limitations  of  available  equipment  precluded  measurement  of  a  complete 
radial  strain  profile  during  high  pressure  testing.  Instead,  the  approach  was  to  derive  a  shape 
function  for  the  bulged  shape  that  could  be  defined  with  a  few  measurements.  The  shape 
function  was  derived  from  a  set  of  passive  radial  deformation  measurements  on  triaxial  test 
specimens.  For  the  passive  measurements,  lines  were  scribed  around  the  circumference  of  a 
virgin  cylindrical  sample  of  Salem  limestone  at  approximately  6  mm  spacings.  The  initial 
diameter  and  height  above  the  base  of  each  line  were  recorded  prior  to  testing.  The  sample  was 
then  subjected  to  a  typical  triaxial  compression  test,  and  the  circumfrential  lines  were 
remeasured.  Active  instrumentation  consisted  of  three  radial  measurements,  one  at  mid-height 
and  two  at  different  heights  on  one  side  of  the  middle,  instead  of  the  usual  single  gage  at  mid¬ 
height. 


The  passive  radial  deformation  data  obtained  over  the  full  sample  height  were  used  to 
derive  a  shape  function  which,  in  turn,  was  used  to  derive  total  volume  strain  estimates  at  each 
time  increment  from  the  three  radial  deformation  measurements.  As  a  validation  of  the  volume 
strain  measurement  procedure,  comparisons  were  made  against  total  volume  strain  estimates 
derived  from  pore  fluid  displacement  measurements  made  during  saturated  drained  tests.  Here, 
the  volume  of  water  drained  from  the  sample  during  testing  was  measured  and  used  to  compute 
an  alternative  estimate  of  total  volume  strain.  Throughout  this  procedure,  the  volume 
determination  based  on  the  shape  function  and  three  radial  displacement  records  was  found  to 
be  imperfect,  but  a  very  great  improvement  over  the  conventional  technique. 

Stress-strain  data  from  a  typical  200  MPa  undrained  triaxial  compression  test  are  shown 
in  Figure  5.21.  The  figure  does  not  include  the  hydrostatic  loading  portion  of  the  test  and 
begins  with  the  initiation  of  the  shear  phase,  i.e.,  the  application  of  the  axial  strain  at  a  constant 
confining  pressure.  The  axial  strain  plotted  in  Figure  5.21  is  an  average  over  the  entire  sample 
height.  Three  radial  strain  curves  are  presented,  corresponding  to  three  different  positions  as 
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shown  in  the  figure.  The  stress  strain  response  is  approximately  linear  up  to  an  axial  strain 
value  of  0.95%  where  brittle  failure  occurs  at  a  true  stress  difference  of  125  MPa.  As  axial 
strain  is  further  increased,  there  is  a  slight  drop  in  stress  difference,  followed  by  a  region  of 
approximately  constant  stress.  The  stress  difference  then  increases  slightly  at  strains  beyond 
5%,  eventually  reaching  almost  140  MPa  before  the  test  was  unloaded  at  15%  axial  strain. 

Figure  5.22  depicts  the  measured  pore  pressure  and  Figure  5.23  presents  volume  strains 
determined  through  two  different  techniques  from  the  same  200  MPa  undrained  test  shown  in 
Figure  5.21.  Axial  strain  was  selected  as  the  independent  variable  for  these  |}'ots  because  it  is 
independently  ccntrolled  during  testing.  As  shown  in  Figure  5.22,  the  pore  picssure  increases 
from  105  MPa  at  the  end  of  the  hydrostatic  phase  to  145  MPa  at  the  point  of  brittle  failure 
(0.95%  axial  strain).  From  that  point,  there  is  a  slight  additional  increase  in  pore  pressure 
followed  by  a  very  slight  decrease  to  142  MPa  at  15%  axial  strain.  During  the  post  failure 
shear,  the  pore  pressure  is  nearly  constant. 

The  volume  strain  determined  in  the  conventional  manner  using  the  radial  strain  measured 
at  mid-height  and  the  average  axial  strain.  Figure  5.23,  shows  an  initial  compressive  volume 
strain  peaking  when  the  axial  strain  reaches  0.95%,  corresponding  to  brittle  failure  of  the 
sample.  As  axial  strain  increases,  that  measure  of  volume  strain  decreases,  i.e.,  becomes  more 
dilatant,  until  it  reaches  a  value  of  approximately  -7.5%  at  15%  axial  strain.  Comparison  in 
Figure  5.20  of  the  relationship  between  pore  pressure  and  volume  strain  before  and  after  the 
brittle  failure  point  leads  to  the  conclusion  that  the  apparent  large  increase  in  volume  suggested 
by  the  conventional  volume  strain  determination  is  inconsistent  with  the  measured  pore  pressure. 
This  result  is  to  be  expected  since  the  conventional  volume  strain  is  derived  from  a  radial 
deformation  measurement  at  the  largest  point  in  a  nonuniform  radial  strain  distribution,  a 
procedure  which  clearly  overestimates  the  volume  strain  relative  to  the  average  total  volume 
strain. 


The  total  volume  strain  determination,  based  on  three  radial  deformation  measurements 
and  the  shape  function,  the  second  curve  presented  in  Figure  5.23  shows  markedly  different 
behavior.  In  the  small  strain  region,  prior  to  brittle  failure,  the  two  curves  are  essentially 
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identical.  Beyond  that  point,  the  total  volume  strain  line  stays  within  0.5%  of  the  zero  volume 
strain  line.  While  this  is  much  more  consistent  with  the  pore  pressure  behavior,  the  variations 
are  still  greater  than  expected.  Note  that  the  variation  in  total  volumetric  strain  is  greater  after 
failure  (0.95%  axial  strain)  than  before  but  the  pore  pressure  change  occurs  almost  entirely 
before.  It  appears  that  both  measures  of  volume  strain  result  in  a  reasonably  correct 
approximation  in  the  region  prior  to  brittle  failure  where  the  strains  are  relatively  small  and 
uniform.  Subsequent  to  failure,  both  axial  and  radial  strains  grow  very  large  with  essentially 
zero  change  in  stress  difference.  Here,  it  is  very  difficult  to  obtain  an  accurate  measure  of 
volume  strain  and  neither  method  gives  a  volume  strain  determination  with  less  than  0.5% 
variation  in  volume  strain  that  would  be  required  to  completely  resolve  the  details  of  pore 
pressure  and  volume  strain  behavior. 

In  the  course  of  performing  the  passive  shape  measurements  on  the  triaxial  samples,  it 
became  apparent  that  there  are  nonuniform  axial,  as  well  as  radial,  deformations.  Based  on  the 
pre-  and  post-test  measurements,  axial  and  radial  strains  were  computed  for  each  (approximately 
6  mm)  segment  of  the  sample.  The  passive  measurements  were  not  designed  to  produce 
accurate  axial  deformation  information.  The  width  of  the  scribed  lines  was  approximately  0.5 
mm,  and  over  the  gage  length  of  6  mm,  an  error  of  half  a  line  width  translates  to  5%  axial 
strain.  In  spite  of  the  measurement  inaccuracies,  the  passive  data  clearly  indicate  the  axial  strain 
is  larger  at  the  center  of  the  sample  than  at  the  ends.  Figure  5.24  shows  profiles  of  radial,  axial 
and  volume  strain  along  the  axis  of  a  47.5  mm  diameter  by  95.8  mm  long  sample  that  were 
computed  from  the  passive  measurements,  along  with  an  indication  of  the  average  axial  strain 
based  on  the  overall  length  change  of  the  sample.  Figure  5.25  presents  a  scatter  plot  of  radial 
and  axial  strain  computed  from  the  passive  measurements  for  the  individual  sample  segments. 
Also  shown  is  a  least  squares  fit  demonstrating  the  trend  in  the  data. 

The  non-uniformity  of  axial  strain  is  significant  even  in  the  case  of  drained  triaxial 
compression  where  the  conventional  assumptions  described  in  Section  5.3  are  applicable.  To 
the  first  order,  the  volume  strain  computation  consists  of  adding  the  axial  strain,  which  is 
positive,  and  twice  the  radial  strain,  which  is  negative.  Thus,  the  volume  strain  is  the  difference 
of  two  numbers  of  like  magnitude  and  is  very  sensitive  to  small  changes  in  either  number.  In 
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the  case  where  the  test  section  is  considered  to  be  the  central  third  of  the  sample  and  the  radial 
strain  is  measured  at  or  near  the  largest  point  on  the  bulged  section,  an  underestimate  of  the 
axial  strain  near  the  mid  height  would  result  in  an  overestimate  of  the  dilatancy  at  that 
measurement  point.  This  is  exactly  what  happens  when  the  axial  strain  is  measured  in  the 
conventional  manner  by  averaging  over  the  entire  length  of  the  specimen. 


5.4  NUMERICAL  SIMULATIONS  OF  UIVDRAINED  TRIAXIAL  COMPRESSION 

TESTS 

As  described  in  the  previous  section,  it  is  known  that  the  loading  conditions  in  real 
triaxial  compression  tests,  while  approximating  the  desired  perfectly  normal  traction,  actually 
include  shear  stresses  that  cause  significant  non-uniformities  in  the  induced  stress  and  strain 
fields.  In  an  effort  to  better  understand  the  actual  behavior  of  the  material  during  these  tests, 
finite  element  simulations  of  the  test  were  performed  using  MEM,  the  two-phase  material 
modeling  program  with  fluid  transport  that  is  described  in  Section  2.  The  material  models 
discussed  in  Section  3  and  the  parameters  discussed  earlier  in  this  section  were  used  in  these 
calculations.  An  effort  was  made  to  correctly  model  the  non-ideal  behavior  that  occurs  during 
a  real  test.  A  ten-by-ten  multi-element  axisymmetric  grid  was  used  in  order  to  model  any  stress 
and  strain  gradients  in  the  sample.  A  symmetry  boundary  condition  was  used  at  sample  mid¬ 
height,  making  it  necessary  to  explicitly  model  only  one  half  of  the  sample.  The  finite  element 
grid  and  boundary  conditions  are  shown  in  Figure  5.26.  At  the  end  where  the  hardened  steel 
end  cap  bears  on  the  sample,  boundary  constraints  were  set  to  allow  no  deformation  parallel  to 
the  interface,  an  idealization  that  slightly  exaggerates  the  physical  effect  observed  in  the  tests. 
The  flow  of  pore  water  from  element  to  element  within  the  mesh  was  modeled  in  these 
simulations. 

Since  the  multi-element  version  of  MEM  does  not  currently  have  the  capability  to  change 
from  stress  to  displacement  boundary  conditions  in  the  course  of  one  run,  the  analysis  was 
started  with  the  hydrostatic  load  already  in  place.  The  initial  pore  pressure,  porosity,  and 
permeability  were  computed  based  on  the  200  MPa  confining  pressure  and  input  to  the 
calculation.  The  confining  pressure  was  held  constant  throughout  the  calculation  and  no  pore 
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fluid  drainage  was  allowed.  In  the  axial  direction  the  loading  was  strain  controlled.  The  applied 
axial  strain  was  linearly  increased  with  time  to  a  maximum  of  15%  at  3000  seconds  ^50 
minutes),  approximating  the  duration  of  an  actual  laboratory  test. 

Figure  5.27  presents  radial,  axial,  and  volume  strain  contours  computed  from  the  output 
of  the  MEM  analysis  of  a  200  MPa  undrained  triaxial  compression  test.  The  plots  in  Figure 
5.27  represent  the  end  of  the  analysis  where  the  average  axial  strain  was  15%  (positive  strain 
represents  compression).  The  contour  plots  exhibit  the  same  trends  observed  in  the  test  data. 
There  is  a  significant  variation  in  both  axial  and  radial  strain  along  the  axis  of  the  sample.  The 
axial  strain  varies  from  4  to  10%  at  the  end  of  the  sample  and  from  20  to  26%  at  mid-height. 
At  approximately  the  quarter  point  of  the  sample,  the  axial  strain  is  equal  to  the  average  value 
of  15%.  Similarly,  in  the  radial  direction  the  strain  is  essentially  zero  at  the  end  cap,  where  the 
boundary  condition  prevents  radial  deformation,  and  at  mid-height  radial  expansion  ranges  from 
10  to  15%.  Taken  together  to  arrive  at  volume  strains,  these  data  show  a  net  compression  near 
the  end  of  the  sample  and  a  net  expansion  at  the  center.  Another  result  that  has  been 
consistently  obtained  in  the  calculations,  but  which  does  not  make  a  very  interesting  plot  is 
constant  pore  pressure  throughout  the  grid.  This  indicates  that  strain  rates  during  the  test  are 
slow  enough  to  permit  pore  pressures  to  equilibrate  and  that  pore  pressures  measured  at  the  base 
are  representative  of  those  throughout  the  specimen. 

The  above  results  are  all  at  least  qualitively  in  agreement  with  the  test  data.  The  volume 
strain  contours  clearly  show  that  the  sample  undergoes  a  loss  of  void  volume  near  the  end  caps 
with  a  corresponding  void  volume  increase  near  the  center.  This  is  accompanied  by  a 
redistribution  of  pore  water  within  the  sample,  with  the  '-  ater  migrating  from  the  ends  toward 
the  center. 

One  way  of  comparing  the  MEM  results  with  test  data  is  to  assume  the  calculational 
results  represent  test  measurement  and  to  process  them  in  exactly  the  same  manner  as  test  data. 
Figure  5.28  shows  pore  pressure  plotted  as  a  function  of  volume  strain.  Included  are  two 
different  curves  computed  from  the  200  MPa  MEM  simulation  along  with  the  same  test  data  that 
was  presented  in  Figure  5.20.  The  dashed  curve  was  computed  from  the  MEM  output  in  the 


s^me  way  the  test  data  curve  was  computed  from  measurements.  In  boih  cases  the  conventional 
approach  was  used,  with  axial  strain  based  on  overall  axial  deformation  of  the  sample  and  radial 
strain  computed  from  the  radial  deformation  at  mid-height.  Both  curves  indicate  that  the  pore 
pressure  init-ally  increases  to  a  peak,  corresponding  to  the  initial  compression  and  cementation 
crushing,  at  a  fraction  of  a  percent  positive  (compressive)  volume  strain.  In  the  initial  loading 
regime  where  volume  strains  are  relatively  uniform  and  measurements  are  relative' y  accurate, 
the  slopes  of  the  curves  representing  test  data  and  finite  element  analysis  are  reasonably 
consistent  and  are  representative  of  the  expected  sensitivity  of  pore  oressure  to  volume  strain. 

As  loading  continues,  the  pore  pressure  in  the  test  data  remains  essentially  constant  while 
there  is  a  steady  pore  pressure  decrease  in  the  calculation.  This  indicates  that  the  net  pt.c 
volume  in  the  test  sample  stayed  approximately  constant  through  this  phase  of  loading,  ’n 
contrast,  the  calculation  appears  to  have  undergone  a  slight  volumetric  expansion,  resulting  in 
a  pore  pressure  drop.  The  dilatancy  in  the  computational  model  of  the  rock  skeleton  is  not 
surprising  considering  the  fact  that  the  skeleton  model  was  derived  from  test  data  processed  in 
the  conventional  manner  which  has  been  shown  to  overestimate  dilatancy. 

Both  the  test  data  and  the  MEM  output  that  was  processed  like  test  data  have  rates  of 
change  of  pore  pressure  with  volume  strain  in  the  post-crush-up  regime  that  are  significantly 
low,  i.e.,  well  below  the  expected  sensitivity  of  pore  pressure  to  volume  strain  as  characterized 
by  comparison  with  the  initial  loading  slope.  The  third  curve  on  Figure  5.28  shows  volume 
strain  computed  from  the  MEM  output  by  taking  the  actual  volume  change  element  by  element 
and  summing  over  the  entire  grid.  The  slope  of  this  curve  is  consistent  with  the  initial  loading 
slope  as  the  skeleton  dilates  and  the  pore  pressure  diminishes. 

Figure  5.29  is  a  plot  of  volume  strain  as  a  function  of  axial  strain,  comparing 
experimental  and  analytical  results  for  both  the  volume  strain  determined  by  the  conventional 
laboratory  approach  and  the  total  volume  strain  averaged  over  the  enMre  sample.  The  figure 
shows  that  when  the  volume  strains  are  computed  consistantly,  there  is  ^ood  agreement  between 
the  test  data  and  the  finite  element  simulation.  The  figure  also  clearly  demonstrates  the 
importance  of  correctly  interpreting  the  strain  definitions  in  a  manner  relevant  to  undrained 
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testing.  Finally,  comparison  of  Figure  5.29  to  5.28  emphasizes  the  sensitivity  of  the  pore 
pressure  to  changes  in  volume  strain  and  suggests  that  pore  pressure  response  data  from 
undrained  tests  might  provide  an  excellent  means  of  fine  tuning  the  dilatancy  parameters  in  the 
material  model. 


5.5  CONCLUSIONS  RELATED  TO  UNDRAINED  TRIAXIAL  TESTING 

The  use  of  numerical  calculations  of  laboratory  tests  is  a  powerful  tool  for  improved 
understanding  of  test  results  and  for  obtaining  more  accurate  material  model  parameters.  The 
comparisons  presented  in  this  section  have  demonstrated  good  agreement  between  laboratory  test 
results  and  two-phase  finite  element  simulation  with  fluid  transport.  Because  of  inaccuracies  in 
the  skeleton  model  parameters,  the  simulation  predicted  more  dilatancy  than  was  actually 
observed  in  the  laboratory.  This  raises  the  possibility  of  performing  iterative  numerical 
calculations  of  test  results  to  better  define  the  model  parameters.  Since  the  pore  pressure 
response  in  an  undrained  test  is  very  sensitive  to  small  changes  in  pore  volume,  pore  pressure 
response  is  well  suited  for  fine  tuning  the  volumetric  strain  behavior  of  a  model. 

The  lab  test  results  and  the  simulations  are  both  consistent  with  the  fundamental  view  that 
pore  pressure  response  is  a  function  of  total  volume  change.  In  an  undrained  triaxial 
compression  test,  volumetric  compression  occurs  near  the  ends  of  the  sample  because  of  the 
lateral  constraint  imposed  by  the  end  caps.  Near  the  center  of  the  sample,  volumetric  expansion 
occurs  due  to  shear  induced  dilatancy.  Because  of  the  variation  in  volumetric  response  along 
the  axis  of  the  sample,  pore  pressure  gradients  are  induced  which  in  turn  cause  flow  of  pore 
water  from  the  ends  of  the  sample  where  the  volume  is  being  reduced  to  the  center  of  the  sample 
which  is  expanding. 

While  it  is  clear  from  post-test  inspection  of  triaxial  compression  test  samples  that  there 
is  a  non-uniform  radial  strain  response,  both  laboratory  measurements  and  finite  element 
simulations  support  the  less  obvious  conclusion  that  there  are  also  non-uniform  axial  strains. 
Axial  strains  are  smaller  than  average  in  the  constrained  material  near  the  ends  of  the  sample 
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and  larger  than  average  in  the  less  constrained  material  at  the  center  of  the  sample.  Thus,  the 
volume  change  response  computed  at  the  center  of  the  sample  based  on  radial  strain  at  the  center 
and  average  axial  strain  tends  to  be  significantly  more  expansive  than  the  actual  behavior  at  any 
point  in  the  sample.  This  introduces  errors  in  model  parameters  causing  the  models  to  predict 
more  dilatancy  than  actually  occurs. 
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Table  5.1.  Physical  Properties  of  Salem  limestone. 


Property 


Value 


Units 


Specific  Gravity  of  the  Solid  Grains 


G, 


2.72 


Porosity 


0.128 


Saturation 


100 


% 


Dry  Bulk  Density 


Tdr 


Saturated  Bulk  Density 


"Ysat 


2.37 


2.50 


Mg/m^ 


Mg/m^ 


Table  5.2.  Grain  compressibility  parameters  for  Salem  limestone. 


Parameter 

Value 

Units 

Initial  bulk  modulus  of  grains 

Ks 

69,000 

MPa 

Initial  mass  density  of  grains 

Po 

2.72 

Mg/m^ 

Initial  wave  velocity  at  low  pressure 

Co 

6,760 

m/s 

Initial  Poisson’s  ratio 

0.25 

— 

Constant  relating  loading  wave  velocity  to  peak  particle  velocity 

s 

1.5 

— 

Threshold  pressure  beyond  which  solid  grains  behave  like  a  fluid 

Pb 

5,000 

MPa 
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Figure  5.1.  Pressure  vs.  volumetric  strain  in  cycled  hydrostatic  com])ression  test  for 
Salem  limestone. 
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Axial  Strain  (%) 

Figure  5.2b.  Volumetric  strain  vs.  axial  strain  in  triaxial  compression  test  for  Salem  Limestone. 
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Figure  5.3.  Elastic  Young's  modulus  as  a  function  of  previous  maximum  mean  stress 
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Figure  5.4.  Compressive  plastic  work  as  a  function  of  mean  stress  In  the  virgin  liydrostatic 
compression  test. 
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Figure  5.5.  Shear  screngtli  enveJojie. 
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Figure  5.6. 


HOZ  lit. 
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Expansive  plastic  hardening  function  vs.  normalized  expansive  plastic  work 
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Figure  5.12.  Coefficient  t.  as  a  function  of  mean  stress 
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Figure  5.15b.  Volumetric  strain  vs.  axial  strain  in  the  drained  trlaxial  compression  at 
On  =  50  MP  ,  MEM  prediction  compared  with  lal)  measurement. 
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Figure  5.16b. 
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Figure  5.17b.  Vol 
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Figure  5.19.  Axial  stress  vs.  axial  strain  in  the  drained  uniaxi.il  strain  loading 
compared  with  lab  measurement. 
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Figure  5.22. 
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Figure  5.23.  Total  volume  strain  determined  from  three  active  radial  measurements  and  the  sampl 
shape  function  compared  to  conventional  volume  strain  (2  kb  undrained  triaxial 
compression  of  Salem  limestone). 
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Figure  5.28.  Comparison  of  pore  pressure  response  computed  using  MKM  wi.'.i  test  data  for  trlaxial 
loading  of  Salem  limestone  at  2  kb  confining  pressure. 


SECTION  6 

DEVELOPMENT  OF  EQUIPMENT  TO  PERFORM  WAVE  PROPAGATION 
STUDIES  IN  SATURATED  POROUS  MEDIA 


Both  theoretical  and  analytical  solutions  for  wave  propagation  in  saturated  porous 
media  were  presented  by  Kim,  Blouin,  and  Timian  (1987),  and  expanded  by  Kim,  Blouin, 
Chitty,  and  Merkle  (1988).  These  solutions  predict  two  types  of  waves,  termed  waves  of  the 
first  and  second  kinds.  The  wave  of  the  first  kind  appears  to  correspond  to  the  familiar 
compressional  wave  in  a  single  phase  medium.  It  is  characterized  by  the  skeleton  and  pore 
fluid,  both  in  compression,  moving  in  the  direction  of  propagation.  The  wave  of  the  second 
kind  propagates  at  a  lower  velocity  with  compression  in  the  pore  fluid  and  tensile  effective 
stress  in  the  skeleton.  These  analyses  also  predict  that  certain  aspects  of  the  wave 
propagation  behavior  are  dependent  on  the  frequency  of  excitation,  i.e.  pulse  duration,  as 
illustrated  in  Figure  6.1.  The  speed  of  the  wave  of  the  first  kind  varies  with  excitation 
frequency,  asymptotically  approaching  a  lower  bound  at  low  frequencies,  an  upper  bound  at 
high  frequencies  and  a  having  transition  in  the  intermediate  frequency  range.  For  the 
particular  combination  of  materials  considered  in  the  analysis,  the  transition  occurs  between 
about  10,000  and  100,000  rad/s.  The  energy  dissipation  in  the  wave  of  the  first  kind  is 
significantly  greater  in  the  transition  region  than  at  frequencies  either  above  or  below. 
Similarly,  the  wavespeed  and  energy  dissipation  of  the  wave  of  the  second  kind  are  dependent 
on  excitation  frequency  as  illustrated  in  Figure  6.1. 

While  the  analytical  and  numerical  solutions  consistently  predict  these  phenomena, 
some  of  them  have  not  been  observed  in  the  laboratory.  Thus,  in  order  to  verify  the  analysis 
results,  and  gain  additional  insight  into  the  physical  processes  involved,  a  program  of 
laboratory  experiments  was  undertaken.  This  involved  construction  of  a  laboratory  apparatus 
capable  of  propagating  compressive  waves  through  saturated  porous  materials  under 
controlled  conditions.  This  section  of  the  report  describes  the  design  of  laboratory  test 
apparatus  and  presents  the  results  of  a  series  of  parameter  calculations  that  were  used  to  guide 
the  equipment  design. 
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The  calculations  presented  in  this  section  of  the  rt^rt  were  performed  prior  to  the 
availability  of  property  measurements  for  the  porous  stainless  steel  bars  that  were  eventually 
used  as  test  specimens.  The  material  properties  used  in  these  calculations  were  estimated 
based  on  available  product  information  supplied  by  manufacturers  that  turned  out  to  be 
significantly  different  than  the  properties  that  were  eventually  measured.  Section  7  presents 
the  properties  of  the  specimen  that  was  tested  and  the  results  of  analyses  based  on  those 
updated  properties. 

6.1  LABORATORY  TEST  APPARATUS 

The  principal  design  objective  of  the  test  apparatus  was  to  provide  a  means  for 
detailed  study  of  waves  propagating  through  a  fully  saturated  porous  medium.  Specific 
design  objectives  included  requirements  to:  saturate  a  specimen  of  porous  material  with 
controlled  boundary  conditions,  induce  a  compressive  wave  in  the  saturated  specimen  and 
measure  the  resulting  disturbance,  preferably  in  both  the  fluid  and  solid  phases.  The  design 
that  was  selected  for  the  Wave  Propagation  Bar  (WPB)  is  similar  in  appearance  to  a  Split 
Hopkinson  Pressure  Bar  (SHPB),  which  is  described  by  Nichols  (1982).  However,  there  are 
important  differences  between  the  two  devices.  The  SHPB  uses  a  sample  much  shorter  than 
the  wavelength  of  excitation  so  that  there  are  several  reflections  through  the  sample  in  the 
duration  of  the  loading  and  the  strain  field  in  the  sample  is  assumed  to  be  uniform.  Strain 
measurements  are  made  on  metal  bars  on  either  side  of  the  sample  and  a  mathematical 
transformation  is  used  to  infer  the  strain  in  the  sample.  In  contrast,  the  wave  propagation  bar 
being  developed  under  this  project  uses  a  sample  longer  than  the  excitation  wavelength  so 
that  wave  propagation  along  the  sample  can  be  observed  by  strain  measurements  made 
directly  on  the  sample.  Details  of  the  analyses  that  lead  to  :he  specific  design  decisions  were 
presented  by  Blouin,  et  al.  (1990).  The  resulting  design  is  depicted  in  Figure  6.2  and 
described  briefly  in  the  following  paragraphs. 

The  device  is  designed  to  accommodate  a  specimen  of  porous  material  in  the  shape  of 
a  cylinder  2  inches  in  diameter  and  up  to  24  inches  long.  The  specimen  is  held  between  two 
50-inch  long  hardened  steel  bars,  also  of  2-inch  diameter,  designated  the  incident  bar  and 
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transmitter  bar.  On  the  incident  bar  end,  is  a  third  2-inch  diameter  bar,  called  the  striker  bar. 
In  order  to  induce  compressive  waves  in  the  test  specimen,  a  gas  gun  is  used  to  propel  the 
striker  bar  through  a  guide  tube  to  a  point  where  it  impacts  the  incident  bar.  The  impact 
results  in  a  compressive  wave  that  propagates  through  the  incident  bar  and  into  the  specimen. 
The  system  has  been  used  with  striker  bars  ranging  from  5  inches  to  20  inches  long.  This 
corresponds  to  a  range  in  pulse  durations  of  50  to  200  |is.  If  the  pulse  duration  is  assured  to 
be  half  the  period  of  a  sine  wave,  the  resulting  range  of  frequencies  is  2,500  to  10,000  Hz  or 
15,700  to  62,800  rad/s.  This  range  could  be  extended  somewhat  on  either  end  by  fabrication 
of  different  striker  bars. 

The  duration  of  the  excitation  pulse  was  identified  in  the  previous  analytical  work  as  a 
key  variable  responsible  for  changes  in  the  wave  propagation  behavior  that  are  under 
investigation.  In  the  test  apparatus,  the  duration  of  the  incident  wave  is  determined  by  the 
length  of  the  striker  bar.  When  the  striker  bar  impacts  the  incident  bar,  the  time  before  it 
rebounds,  and  hence  the  duration  of  the  incident  pulse,  is  equal  to  the  time  required  for  a 
compressive  wave  to  travel  from  the  point  of  impact  to  the  opposite  end  of  the  striker  bar  and 
return  to  the  impact  point  as  a  tensile  wave  which  causes  separation  between  the  two  bars. 

The  magnitude  of  the  incident  wave  is  determined  by  the  velocity  of  impact,  which  is 
determined  by  the  pressure  in  the  gas  gun.  For  this  work,  the  tests  were  designed  with  the 
objective  of  keeping  all  material  response  in  the  elastic  regime.  Thus,  the  gas  gun  pressure 
was  selected  to  give  adequate  resolution  in  the  instrumentation  without  yielding  the  test 
specimen.  Since  the  material  stays  elastic,  the  exact  magnitude  of  the  stress  wave  does  not 
strongly  influence  the  outcome  of  the  experiments. 

In  order  to  achieve  saturation  of  the  specimen,  all  of  the  air  in  its  pore  space  must  be 
removed  and  replaced  with  fluid.  Further,  in  order  to  apply  the  necessary  boundary 
conditions,  we  impose  the  requirement  that  no  fluid  flow  into  or  out  of  the  specimen  during 
testing.  This  condition  is  accomplished  by  encasing  the  specimen  in  a  Jacket,  a  tight-fitting 
flexible  membrane,  which  is  sealed  to  the  solid  bars  on  either  end  of  the  test  specimen.  The 
exterior  of  the  jacket  is  pressurized  to  a  level  greater  than  the  maximum  pore  pressure, 
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ensuring  that  the  jacket  remains  tight  against  the  test  specimen  and  that  no  fluid  flows  across 
the  specimen  boundary. 

Since  the  exterior  of  the  jacketed  specimen  must  be  pressurized  to  maintain  the 
required  boundary  conditions,  it  is  necessary  that  the  test  apparatus  include  a  pressure  vessel 
to  contain  the  specimen.  In  this  case,  the  pressure  vessel  consists  of  a  standard  double  ended 
hydraulic  cylinder  in  which  the  piston  has  been  replaced  by  the  test  specimen,  with  the 
incident  and  transmitter  bars  passing  through  the  seals  on  either  end  of  the  cylinder. 

In  light  of  the  caiculational  results  presented  in  Section  6.2,  a  steel  sleeve  was  added 
to  the  pressure  vessel,  leaving  a  gap  of  approximately  0.2  inch  between  the  jacket  covering 
the  specimen  and  the  inner  surface  of  the  sleeve.  This  results  in  a  stiffer  radial  confinement 
condition  which  is  intermediate  in  stiffness  between  a  constant-pressure  fluid  confinement  and 
a  uniaxial  strain  condition.  The  latter  is  desirable,  but  is  not  readily  achievable  in  the 
laboratory. 

In  addition  to  the  radial  confinement  provided  by  the  sleeve  and  fluid  in  the  pressure 
vessel,  axial  confinement  is  required  to  maintain  the  required  boundary  conditions.  This  is 
provided  by  a  hydraulic  cylinder  attached  to  the  end  of  the  transmitter  bar.  To  react  against 
the  axial  force  applied  by  the  hydraulic  cylinder,  a  collar  attached  to  the  incident  bar  bears 
against  a  fixed  stop. 

Connection  to  the  pore  space  of  the  installed  test  specimen  is  required  to  evacuate  the 
air  and  introduce  pore  fluid.  This  is  achieved  by  means  of  an  axial  hole  in  the  transmitter  bar 
originating  at  the  specimen  and  intersecting  a  radial  hole  outside  the  pressure  vessel.  The 
entire  apparatus  is  assembled  on  top  of  a  27-inch  deep  wide  flange  steel  beam  which  provides 
a  stiff  surface  to  align  the  parts  and  a  convenient  working  elevation. 
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6J  CONSTRUCTION  OF  THE  TEST  EQUIPMENT 


The  device  described  above  was  constructed  in  ARA’s  Materials  Testing  Laboratory 
in  South  Royalton,  Vermont.  A  20-ft  section  of  27-inch  deep  by  10-inch  wide  steel  I-beam 
(W27X84)  was  drilled  for  attachment  of  the  various  components,  shimmed  level  and  fastened 
to  the  floor  in  the  laboratory.  Both  the  pressure  vessel  and  the  axial  loading  cylinder  were 
fabricated  to  specification  by  a  hydraulic  cylinder  vendor.  The  pressure  vessel  was  retrofitted 
with  electrical  penetrations  to  pass  instrumentation  signals  through  the  pressure  boundary. 

The  pressure  vessel  and  axial  loading  cylinder  were  carefully  aligned  and  bolted  to  the  top  of 
the  steel  beam.  Also  attached  to  the  top  of  the  beam  are  the  guide  tube  for  the  striker  bar 
and  a  bracket  to  react  to  the  axial  load  in  the  bar/specimen  assembly. 

The  gas  gun  incorporated  into  this  apparatus  was  available  at  ARA  at  the  time  of  the 
design.  While  it  is  much  larger  than  required  for  this  purpose,  it  was  the  most  economical 
choice  since  no  additional  hardware  purchase  was  required.  In  order  to  isolate  any 
mechanical  noise  developed  in  the  gas  gun  from  the  instrumented  portion  of  the  test 
assembly,  the  gas  gun  is  supported  on  a  separate  steel  base  structure.  The  support  structure 
for  the  gas  gun  and  the  steel  beam  supporting  the  bar  assembly  are  aligned  so  that  the  outlet 
of  the  gas  gun  slides  a  short  distance  into  the  upstream  end  of  the  striker  bar  guide  tube 
without  any  mechanical  contact. 

Figures  6.3  through  6.7  are  photographs  of  the  completed  test  device.  Figure  6.3 
shows  as  much  as  possible  of  the  entire  assembly,  which  is  over  20  ft  long.  The  gas  gun  is 
out  of  the  picture  to  the  right.  On  the  right  end  of  the  picture  is  the  guide  tube  for  the  striker 
bar.  Moving  to  the  left  are  the  incident  bar,  pressure  vessel,  transmitter  bar,  and  the  axial 
loading  cylinder.  Figure  6.4  presents  a  different  view  of  the  assembly,  looking  down  axis 
with  guide  tube  in  the  foreground.  The  upstream  end  of  the  incident  bar  is  shown  in  Figure 
6.5.  Also  shown  in  Figure  6.5  are  the  downstream  end  of  the  striker  bar  guide  tube,  the 
collar  and  reaction  frame  that  support  the  axial  load  in  the  incident  bar,  and'  a  strain  gage 
station  near  the  upstream  end  of  the  incident  bar.  Figure  6.6  shows  the  incident  bar  entering 
the  pressure  vessel  that  contains  the  test  specimen.  The  hose  and  valves  are  for  filling  and 
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draining  the  pressure  vessel.  The  downstream  end  of  the  pressure  vessel  is  shown  in  Figure 
6.7,  with  the  transmitter  bar  passing  through  a  seal  in  the  end  of  the  vessel.  Also  on  this  end 
of  the  pressure  vessel  are  the  electrical  penetrations  used  to  transmit  signals  from  the 
instrumentation  internal  to  the  pressure  vessel  out  to  the  recording  equipment.  Attached  to 
the  transmitter  bar,  near  the  seal  of  the  pressure  vessel,  is  a  tubing  fitting  with  integral  valve 
which  connects  to  the  pore  space  in  the  test  specimen  by  means  of  an  axial  hole  in  the  end  of 
the  transmitter  bar. 

6.3  PERFORMANCE  TESTS 

A  series  of  performcince  tests  were  conducted  prior  to  testing  a  saturated  porous 
specimen.  For  these  tests  the  device  was  set  up  with  a  2-inch  diameter  by  24-inch  long 
hardwood  specimen  in  place  of  the  planned  saturated  porous  material.  As  with  the  saturated 
porous  metal,  the  hardwood  specimen  has  a  lower  acoustic  impedance  than  the  steel  incident 
bar.  The  tests  were  conducted  with  no  fluid  and  no  confining  pressure  in  the  pressure  vessel. 
The  axial  loading  cylinder  was  pressurized  to  hold  the  incident  bar,  the  hardwood  specimen, 
and  the  transmitter  bar  tightly  together. 

A  pair  of  electrodes  with  one  inch  separation  along  the  axis  of  the  bar  was  installed  at 
the  end  of  the  striker  bar  guide  tube.  Using  a  simple  two-step  voltage  divider  circuit,  and  a 
20  MHz  oscilloscope,  measurements  were  made  of  the  striker  bar  velocity  immediately 
preceding  impact  with  the  incident  bar.  Using  this  procedure,  it  was  determined  that  12  psi 
pressure  in  the  gas  gun  results  in  an  impact  velocity  of  13  -  14  ft/s.  For  like  materials 
impacting,  the  magnitude  of  the  resulting  stress  wave,  a„  is  given  by: 


where:  p  =  mass  density  of  the  material 

Cp  =  compressive  wavespeed  of  the  material 
V  =  impact  velocity 
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A  thin  steel  bar  has  mass  density  and  compiessional  wavespeed  of  0.00073  lb  sVin.  and 
200,000  in./s,  respectively.  At  an  impact  velocity  of  14  ft/s  (=  168  in7s),  the  resulting 
computed  stress  is  12,300  psi.  Similar  velocities  were  obtained  using  3.5  psi  gas  gun 
pressure  with  a  5-inch  striker  bar  and  5  psi  with  a  10-inch  striker  bar. 

Several  trial  runs  were  made  with  12  psi  gas  gun  pressure,  the  20-inch  striker  and  the 
hardwood  specimen.  There  were  no  instruments  on  the  specimen,  but  two  full  strain  gage 
bridges  were  installed  on  the  incident  bar,  one  10.5  inches  from  the  end  that  is  impacted  by 
the  striker  and  the  other  2.25  inches  from  the  specimen  end,  resulting  in  a  separation  distance 
of  37.25  inches.  Both  strain  gage  bridges  were  monitored  with  a  20  MHz  digital  storage 
oscilloscope.  The  data  recorded  in  one  typical  run  are  presented  in  Figure  6.8  where  the 
origin  of  the  time  scale  is  arbitrary.  The  measured  voltages  were  converted  to  stresses  using 
a  calibration  factor  computed  from  the  nominal  values  of  strain  gage  factor,  excitation 
voltage,  amplifier  gain,  elastic  modulus,  and  Poisson’s  ratio.  The  measured  stress  wave 
amplitude  of  approximately  12,300  psi  matches  the  value  computed  from  the  impact  velocity. 

The  duration  of  the  wave  should  be  equal  to  two  wave  transit  times  of  the  striker  bar. 
For  the  20-inch  striker  bar,  the  expected  pulse  duration  is  200  ps.  In  the  measured  signals, 
the  edges  of  the  pulse  exhibit  some  spreading  due  to  dispersion  of  the  wave.  However,  the 
main  part  of  the  wave  very  closely  matches  the  expected  value.  The  pulse  velocity  along  the 
incident  bar  can  be  measured  by  comparing  the  signals  from  the  two  sUain  gage  stations  in 
Figure  6.8.  There  is  more  spreading  on  the  leading  edge  of  the  pulse  at  the  second  strain 
gage  station  than  at  the  first.  If  the  propagation  time  between  the  two  station  is  measured  on 
the  steep  portion  of  the  leading  edge  of  the  waveform,  the  resulting  velocity  is  199,000  in/s 
which  matches  the  value  computed  using  the  (unconflned)  elastic  modulus  of  steel. 
Theoretically,  this  is  the  wave  propagation  velocity  in  a  thin  rod  that  is  free  to  deform 
laterally  as  the  wave  propagates  along  its  axis.  There  is  also  a  component  of  the  wave  that  is 
travelling  faster  as  evidenced  by  the  difference  in  first  arrival  times  at  the  two  stations. 

Based  on  first  arrivals,  a  wavespeed  of  213,000  in/s  is  computed.  This  is  somewhat  ulower 
than  the  231,000  in/s  speed  that  is  computed  for  the  fully  confined  (or  halfspace)  case,  but 
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apparently  indicates  some  tendency  in  that  direction  resulting  from  outrunning  along  the 
central  axis  of  the  bar. 

6.4  PARAMETRIC  CALCULATIONS  OF  LABORATORY  TEST  APPARATUS 

The  propagation  of  stress  waves  in  a  multiphase,  cylindrical  bar  is  a  complex 
mechanical  phenomenon.  In  order  to  design  laboratory  experiments  which  can  be  effectively 
correlated  with  theoretical  predictions,  we  performed  a  series  of  numerical  simulations  of 
various  test  specimen  loading  and  confinement  conditions  to  ensure  that  the  experiments 
would  have  a  maximum  chance  of  success.  These  calculations  were  performed  using  the 
finite  element  code  MPDAP,  described  in  Section  2,  and  the  material  models  discussed  in 
Section  3.  These  calculations  simulate  stress  wave  propagation  in  a  cylindrical  bar  of 
saturated,  porous,  elastic  material  with  the  primary  objective  of  assessing  various  aspects  of 
the  planned  experiments. 

Three  sets  of  parametric  calculations  were  performed  as  shown  in  Table  6.1.  In  the 
first  set  of  calculations  the  porous  bar  is  loaded  in  uniaxial  stress  while  confined  in  a  fluid  at 
constant  confining  pressure.  Both  the  loaded  surface  at  the  end  of  the  bar  and  the  sides  of 
the  bar  are  sealed  to  prevent  pore  fluid  from  moving  across  these  boundaries.  In  the 
experiment  this  would  be  accomplished  by  jacketing  the  sample  between  the  incident  and 
transmitter  bars  within  a  constant  pressure  gas.  Two  types  of  uniaxial  stress  calculations  were 
performed.  In  the  first,  the  pore  fluid  was  free  to  displace  within  the  porous  skeleton  as 
governed  by  a  realistic  permeability  for  the  porous  bar.  In  the  second,  the  movement  of  pore 
fluid  within  the  skeleton  was  prevented  by  using  a  very  low  value  of  permeability.  The 
results  of  the  uniaxial  stress  loading  calculations  showed  very  little  influence  of  relative  pore 
fluid  motion  because,  in  both  the  permeable  and  impermeable  cases,  excess  dynamic  pore 
pressures  within  the  bar  are  relieved  almost  immediately  by  relief  waves  propagating  into  the 
bar  from  the  cylindrical  surface.  Within  the  elastic  regime,  this  phenomenon  was  found  to  be 
independent  of  the  confining  pressure. 
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A  second  set  of  calculations  was  performed  in  which  the  porous  bar  was  loaded  in 
uniaxial  strain.  Based  on  previous  analyses  (Kim,  et  al.,  1988),  it  was  expected  that  this 
loading  condition  would  show  a  dramatic  influence  of  the  relative  pore  fluid  motion.  This 
was  verified  by  comparing  the  permeable  bar  calculation  to  the  companion  calculation  using 
an  impermeable  bar.  Experimentally,  however,  it  is  extremely  difficult  to  achieve  a  true 
dynamic  uniaxial  strain  condition  in  the  laboratory. 

In  order  to  approximate  the  uniaxial  strain  condition  as  closely  as  possible,  so  as  to 
maximize  the  influence  of  relative  pore  fluid  motion,  an  experimental  configuration  using  a 
porous  bar  confined  in  a  small  volume  of  liquid  within  a  rigid  walled  pressure  vessel  was 
considered.  The  idea  of  this  approach  is  to  develop  dynamic  confinement  of  the  porous  bar 
in  the  relatively  stiff  confining  liquid.  The  thin  annulus  of  confining  liquid  provides 
confinement  with  minimum  shear  resistance  and  high  lateral  stiffness.  The  third  set  of 
calculations  replicates  the  experimental  conditions  using  this  test  configuration.  Two 
companion  calculations  were  performed.  The  first  accurately  modeled  the  permeability  of  the 
porous  bar  and  the  second  assumed  that  the  bar  was  impermeable  so  as  to  delineate  the 
influence  of  the  relative  pore  water  motion. 

As  described  by  Blouin,  et  al.  (1990),  porous  stainless  steel  was  selected  for  the  first 
set  of  experiments  because  its  combination  of  high  skeleton  stiffness  and  high  permeability 
should  accentuate  the  influence  of  the  relative  fluid  motion.  The  published  properties  of  the 
materials  chosen  for  the  experiment  were  used  in  these  calculations  and  are  shown  in  Table 
6.2.  The  permeability  properties  of  the  porous  bar  were  obtained  from  the  manufacturer  and 
the  stiffness  properties  were  estimated  based  on  available  information.  Kerosene  was  chosen 
for  the  pore  fluid  because  it  is  non  corrosive. 

Moduli  and  wavespeeds  were  computed  fr^'m  closed  form  solutions  for  the  various 
materials  and  loadings  used  in  the  calculations  widi  the  properties  in  Table  6.2.  These  values 
are  listed  in  Table  6.3.  Both  Young’s  modulus  and  the  constrained  modulus  for  the  porous 
skeleton  are  listed,  as  appropriate,  from  which  wavespeeds  corresponding  to  the  uniaxial 
stress  loading  and  uniaxial  strain  loading,  respectively,  are  computed.  The  wavespeed  for  the 
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uniaxial  stress  loading  of  the  saturated  porous  bar  is  computed  using  Young’s  modulus  of  the 
dry  skeleton  because  of  the  immediate  lateral  relief  of  the  excess  pore  pressures.  Two 
bounding  wavespeeds  are  computed  for  the  uniaxial  strain  loading  of  the  saturated  porous 
steel.  As  described  by  Kim,  et  al.  (1987),  the  compressional  wavespeed  is  a  function  of  the 
permeability  and  the  loading  frequency.  For  low  frequency  excitations  (or  low 
permeabilities),  a  lower  bound  wavespeed  is  obtained  which  is  that  associated  with  the 
undrained  fully  coupled  modulus  computed  according  to  procedures  described  by  Blouin  and 
Kim  (1984).  For  high  frequency  excitations  (or  high  permeabilities),  an  upper  bound 
wavespeed  is  obtained  according  to  the  analytic  procedures  developed  by  Kim,  et  al.  (1987). 
The  computed  wavespeed  relationship  as  a  function  of  frequency  for  the  saturated  porous 
steel  to  be  used  in  the  experimental  validation  is  shown  in  Figure  6.1a.  This  profile  was 
computed  using  the  program  TWAVE  which  is  documented  in  Kim,  et  al.  (1987).  In  part  b 
of  Figure  6.1,  the  damping  for  the  compressional  wave  is  plotted  as  a  function  of  excitation 
frequency.  In  parts  c  and  d,  the  wave  velocity  and  damping  for  waves  of  the  second  kind  are 
plotted  as  functions  of  excitation  frequency. 

A  rectangular  loading  pulse  with  a  peak  stress  of  1,000  psi  was  used  in  all  the 
calculations  listed  in  Table  6.1.  The  loading  pulse  had  a  duration  of  0.1  ms  with  a  rise  time 
and  decay  time  of  I.O  }j.s.  This  loading  excitation  frequency  is  on  the  upper  bound  of  the 
wavespeed  profile  of  Figure  6.1.  As  described  earlier,  the  excitation  frequency  in  the 
laboratory  experiment  can  be  varied  by  changing  the  length  of  the  striker  bar. 

Schematic  section  views  of  the  calculational  grids  and  data  output  locations  for  each 
of  the  three  sets  of  calculations  are  shown  in  Figures  6.9  through  6. 1 1 .  The  uniaxial  stress 
calculations  used  an  axisymmetric  cylindrical  grid  of  1680  elements.  There  are  7  equally 
spaced  elements  in  the  radial  direction  making  up  the  1  inch  radius  porous  bar  and  240 
elements  in  the  axial  direction  making  up  the  24  inch  total  length.  Data  output  locations  are 
shown  at  points  A  through  D  spaced  at  5  inch  increments  within  the  bar.  The  rectangular 
loading  pulse  is  applied  on  the  circular  face  of  the  bar  against  an  impermeable  fiee  surface. 
There  is  an  impermeable  free  boundary  surrounding  the  bar.  In  the  uniaxial  stress 
calculations  reported  here,  no  confining  pressure  was  applied  to  the  lateral  surface,  simulating 
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an  unconfmed  loading.  Other  uniaxial  stress  calculations,  not  reported,  involved  bars 
confined  under  constant  uniform  stress  conditions.  The  response  was  virtually  identical  to 
that  reported  here  for  the  unconfmed  bar. 

The  uniaxial  strain  calculations,  shown  schematically  in  Figure  6.10,  consisted  of  the 
same  rectangular  pressure  pulse  applied  to  a  1  inch  radius  bar  confined  by  an  impermeable 
roller  boundary  such  that  no  lateral  motion  at  the  boundary  is  allowed.  The  element  size  in 
the  axial  direction  is  identical  to  that  used  in  the  uniaxial  stress  calculations. 

The  fluid  confined  calculations,  shown  schematically  in  Figure  6.11,  model  a  1  inch 
radius  saturated  porous  steel  bar  confined  within  a  0.2  inch  annulus  filled  with  kerosene.  The 
kerosene  is  confined  within  a  rigid  impermeable  roller  boundary  which  approximates  the  steel 
pressure  vessel.  The  porous  steel  bar  is  confined  in  the  kerosene  by  an  impermeable  free 
boundary.  The  uniform  rectangular  pressure  loading  is  applied  to  both  the  impermeable 
surface  of  the  bar  and  to  the  confining  fluid.  The  finite  element  grid  for  the  bar  is  identical 
to  that  used  in  the  uniaxial  stress  calculations,  but  an  addidonal  720  elements  are  included  to 
model  the  confining  fluid.  Additional  output  points,  subscripted  with  the  letter  f,  are  included 
to  monitor  the  response  of  the  confining  fluid.  Selected  comparisons  and  analysis  of  the 
calculational  results  are  presented  in  the  following  subsection. 

6,5  ANALYSIS  OF  PARAMETRIC  CALCULATIONS 

/ 

Confinement  conditions  have  a  dramatic  influence  on  the  dynamic  response  of  the 
saturated  porous  bar.  Figure  6.12  compares  the  total  axial  stress  waveforms  at  ranges  of  S, 

10  and  15  inches  from  the  point  of  load  application  (locations  A,  B  and  C  of  Figures  6.9 
through  6.11)  for  the  calculations  in  which  relative  fluid  motion  within  the  porous  skeleton 
was  prevented.  The  waveforms  and  propagation  velocities  of  the  uniaxial  strain  and  fluid 
confined  calculations  are  similar,  but  in  stark  contrast  to  those  of  the  uniaxial  stress 
calculation.  The  rectangular  waveform  of  the  loading  function  is  well  preserved  throughout 
all  ranges  in  the  fluid  confined  and  uniaxial  strain  calculations.  The  wave  in  the  uniaxial 
stress  calculation  is  rapidly  transformed  into  a  sinusoidal  waveform  with  considerable 
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broadening  of  the  pulse  width.  The  propagation  velocity  in  the  uniaxial  strain  calculation,  as 
measured  at  the  half  height  of  the  initial  rise  to  peak  between  the  5  and  15  inch  ranges,  is 
1170  m/s.  This  is  nearly  identical  to  the  calculated  wavespeed  for  the  uniaxial  strain, 
saturated  lower  bound  wavespeed  of  Table  6.3.  As  noted  previously,  this  is  the  applicable 
wavespeed  for  the  case  in  which  there  is  no  relative  fluid  motion.  The  propagation  vek  .  ity 
in  the  uniaxial  stress  calculation,  as  measured  from  the  peak  stress  arrival  between  the  5  and 
15  inch  ranges,  is  835  m/s  which  compares  to  834  m/s  computed  for  the  uniaxial  stress, 
saturated  case  in  Table  6.3. 

The  propagation  velocity  in  the  fluid  confined  porous  bar  is  1210  m/s,  which  is 
slightly  faster  than  that  for  the  uniaxial  strain  loading.  This  indicates  that  the  lateral 
confinement  in  the  fluid  confined  bar  is  somewhat  higher  than  that  in  the  uniaxial  strain 
calculation.  The  propagation  velocities  measured  from  the  three  calculations  with  no  relative 
fluid  motion  are  listed  in  Table  6.4.  Another  indication  of  the  excess  confinement  in  the  fluid 
confined  bar  is  that  the  peak  axial  stresses  are  higher  than  those  in  the  uniaxial  strain 
calculation.  The  total  axial  stresses  are  also  greater  than  the  stress  applied  to  the  end  of  the 
bar. 


Another  interesting  aspect  of  the  propagation  velocity  of  the  uniaxial  stress  pulse  is 
that  the  first  arrival  tends  to  be  coincident  with  the  first  arrival  of  the  uniaxial  strain  pulse. 
This  would  be  expected,  since  some  small  fraction  of  the  energy  propagates  to  the  observation 
point  along  the  center  of  the  bar  at  the  constrained  wavespeed,  uninfluenced  by  the  cylindrical 
unconfined  boundary.  Since  the  bar  is  elastic,  relief  waves  from  the  unconfined  boundary 
cannot  overtake  the  initial  waves  travelling  in  the  interior  of  the  bar. 

The  pore  pressure  profiles,  corresponding  to  the  total  stress  profiles  above,  are  shown 
in  Figure  6.13.  The  shape  of  the  profiles  is  similar  to  that  of  the  total  stresses,  however,  the 
magnitude  of  the  pore  pressure  in  the  uniaxial  stress  case  is  relatively  small,  because  of  the 
relief  provided  at  the  unconfined  boundary  of  the  bar.  The  corresponding  particle  velocity 
and  displacement  profiles  for  the  calculations  with  no  relative  fluid  motion  are  shown  in 
Figures  6.14  and  6.15.  The  velocity  profiles  are  very  similar  in  shape  to  the  total  stress 
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profiles.  The  displacement  profiles  show  larger  displacements  in  the  uniaxial  stress  case,  due 
to  the  lack  of  confinement.  At  early  times,  displacements  in  the  fluid  confined  case  are 
greater  than  those  in  the  uniaxial  case  because  of  the  faster  wavespeed  and  greater  particle 
velocities  in  the  former  case.  However,  at  later  time  the  fluid  confined  and  uniaxial  strain 
displacements  come  into  close  agreement  as  the  significant  negative  velocity  phase  in  the 
fluid  confined  case  works  to  reduce  the  late  time  displacements. 

Comparison  plots  showing  the  influence  of  confinement  for  the  three  sets  of 
calculations  in  which  relative  fluid  motion  is  allowed,  as  governed  by  the  published 
permeability  of  the  porous  steel  bar,  are  shown  in  Figures  6.16  through  6.21.  The  total  stress 
profiles  at  the  5,  10  and  15  inch  ranges  are  shown  in  parts  a,  b  and  c  of  Figure  6.16. 
Comparison  to  the  companion  calculations  with  no  relative  fluid  motion  (Figure  6.12)  shows 
that  there  is  very  little  influence  of  relative  fluid  motion  in  the  uniaxial  stress  case,  but  a 
dramatic  influence  in  the  uniaxial  strain  case.  On  the  basis  of  these  comparisons,  we  decided 
to  replicate  the  uniaxial  strain  condition  as  closely  as  possible  in  our  experiments.  The  total 
stress  profile  for  the  fluid  confined  bar  proposed  for  the  laboratory  experiments  is 
significantly  altered  in  shape  and  reduced  in  amplitude  by  the  relative  fluid  motion,  though 
not  as  dramatically  as  are  the  uniaxial  strain  profiles. 

The  relative  fluid  motion  strongly  influences  the  development  of  early-time  peak  stress 
in  the  uniaxial  strain  loading  as  demonstrated  by  the  comparison  of  the  uniaxial  and  fluid 
confined  waveforms  in  Figure  6.16a.  Only  5  inches  into  the  bar,  the  initially  square  shape  of 
the  wave  has  been  altered  considerably  by  the  relative  fluid  motion.  The  uniaxial  strain 
waveform  is  rapidly  attenuated  and  broadened  compared  to  that  in  the  fluid  confined  bar. 

The  wave  propagation  velocities  measured  from  the  calculations  with  relative  fluid 
motion  are  listed  and  compared  to  the  velocities  without  fluid  motion  in  Table  6.4.  There  is 
little  change  in  the  unconfined  bar  velocity  because  the  relief  from  the  cylindrical  surface  of 
the  bar  tends  to  minimize  pore  pressures  and  pore  pressure  gradients  within  the  bar.  There  is 
a  substantial  increase  in  the  velocity  in  the  bar  under  uniaxial  strain,  from  1170  to  1375  m/s. 
This  is  in  agreement  with  the  velocity  increases  predicted  by  the  analytical  TWAVE 


calculation  of  the  upper  bound  wavespeed  shown  in  Figure  6.1a.  As  explained  by  Kim,  et  al. 
(1987),  the  upper  bound  wavespeed  is  approached  as  either  frequency  or  permeability 
becomes  large.  The  permeability  of  the  bar  should  be  such  that  significantly  higher 
wavespeeds  will  be  measured  than  indicated  by  the  closed  form  solution  for  the  undrained 
fully  coupled  loading  condition  derived  by  Blouin  and  Kim  (1984).  The  wavespeed  in  the 
fluid  confined  bar  also  shows  an  increase  once  relative  motion  is  allowed.  However,  the 
wavespeed  of  1290  m/s  is  85  m/s  slower  than  the  wavespeed  under  pure  uniaxial  strain 
conditions. 

The  pore  pressure  time  history  comparisons  of  Figure  6.17  show  some  trends  similar 
to  those  observed  in  the  total  stress  comparisons,  though  the  initial  rapid  decay  of  the  pore 
pressure  pulse  during  the  uniaxial  strain  loading  is  less  pronounced.  The  pore  pressure  in  the 
uniaxial  strain  calculation  exhibits  more  pulse  broadening  and  modification  than  that  in  the 
fluid  confined  case.  Pore  pressures  developed  in  the  uniaxial  stress  case  are  much  less  than 
those  in  the  other  two  cases  because  of  the  relief  from  the  surface  of  the  bar. 

A  comparison  of  the  skeleton  velocity  time  histories  for  the  three  calculations  with 
relative  fluid  motion  at  various  locations  along  the  bar  is  shown  in  Figure  6.18.  The 
waveforms,  waveform  dispersion,  and  other  alterations  in  shape  closely  match  those  for  the 
stress  time  histories  shown  in  Figure  6.16.  The  peak  velocities  are  significantly  higher  in  the 
uniaxial  stress  loadings  because  of  the  lack  of  confinement.  The  peak  velocity  attenuation  is 
greatest  in  the  uniaxial  strain  case  and  least  in  the  uniaxial  stress  case.  The  corresponding 
displacement  time  histories  are  shown  in  Figure  6.19.  These  reflect  the  differences  in  the 
velocity  time  histories  including  larger  displacement  in  the  uniaxial  stress  case.  There  is 
more  dispersion  in  the  uniaxial  strain  waveform  than  in  that  for  the  fluid  confined  case. 

In  order  to  help  understand  the  differences  in  the  three  calculations  with  relative  fluid 
motion,  as  well  as  the  influence  of  the  relative  motion  on  the  overall  response,  relative  motion 
time  histories  at  three  locations  along  the  axis  of  the  bar  are  presented  in  Figures  6.20  and 
6.21.  Figure  6.20  shows  average  relative  velocity  between  the  pore  fluid  and  the  porous 
skeleton  at  the  5,  10  and  IS  inch  ranges.  The  relative  velocity  waveforms  for  the  uniaxial 
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strain  and  fluid  conflned  cases  are  similar  in  shape,  with  the  fluid  conflned  case  having  the 
highest  amplitude  at  all  ranges.  Waveforms  for  both  conflned  conditions  jump  rapidly  to  a 
peak,  begin  a  gradual  decay,  and  then  drop  suddenly  into  a  negative  phase.  At  the  close  in 
ranges  peak  relative  velocities  are  higher  than  the  corresponding  skeleton  velocities  shown  in 
Figure  6.18.  At  further  ranges  they  are  about  equal  to  the  peak  skeleton  velocities.  Thus, 
just  behind  the  wavefront,  the  pore  fluid  surges  ahead  of  the  skeleton  as  fluid  is  squeezed  out 
of  adjoining  elements  by  compression  of  the  skeleton.  Then  the  relative  velocity  gradually 
diminishes  during  which  time  the  skeleton  velocity  is  gradually  increasing  (see  Figure  6.18). 
This  is  followed  by  a  sharp  drop  in  relative  velocity  and  a  large  negative  phase  during  which 
the  pore  fluid  is  moving  backward  relative  to  the  skeleton.  The  sharp  drop  and  negative 
phase  in  the  relative  fluid  motion  corresponds  with  a  sharp  reduction  in  skeleton  velocity. 

The  end  of  the  negative  phase  occurs  at  about  the  time  the  skeleton  velocity  drops  to  zero. 

The  relative  fluid  velocity  in  the  uniaxial  stress  loading  (Figure  6.20)  is  much  lower  in 
magnitude  and  later  in  time  than  the  corresponding  signal  in  the  confined  bar  calculations. 
This  is  in  keeping  with  previous  analysis  which  showed  that  relief  waves  from  the  sides  of 
the  bar  immediately  relieve  pore  pressures  in  the  fluid.  The  character  of  the  uniaxial  stress 
waveform  is  opposite  of  those  in  the  confined  bar  calculations  in  that  the  relative  velocity 
rises  gradually  to  a  peak  and  then  decays  into  a  negative  phase,  resulting  in  a  sinusoidally 
shaped  waveform. 

The  displacement  of  the  pore  fluid  relative  to  the  porous  skeleton  is  shown  in  Figure 
6.21  for  the  three  ranges  of  interest.  In  the  case  of  the  uniaxial  strain  and  fluid  confined  bar 
calculations  the  pore  fluid  shows  a  strong  surge  ahead  of  the  porous  skeleton  which  peaks 
well  before  the  maximum  skeleton  displacement  (see  Figure  6.19  for  comparison).  By  the 
time  of  peak  skeleton  displacement,  the  pore  fluid  has  moved  back  toward  its  original 
position  in  the  bar.  The  backward  motion  continues  into  a  negative  phase  during  which  the 
pore  fluid  is  behind  its  original  position  in  the  bar.  The  negative  phase  gradually  lessens  until 
the  fluid  returns  to  its  original  position.  There  is  very  little  relative  displacement  in  the 
uniaxial  stress  case,  with  only  a  mild  oscillation  developing. 
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In  the  remaining  illustrations  in  this  section  we  take  a  more  detailed  look  at  the 
influence  of  relative  fluid  motion  on  stress  and  motion  propagation  in  the  uniaxial  stress 
(Figures  6.22  -  6.25)  uniaxial  strain  (Figures  6.26  -  6.29)  and  fluid  confined  loadings  (Figures 
6.30  -  6.33).  Comparison  of  total  stress  time  histories  for  the  uniaxial  stress  loadings  at  the 
three  ranges  of  interest  with  and  without  relative  fluid  motion  are  shown  in  Figure  6.22. 

There  is  very  little  difference  in  the  waveforms  because  so  little  relative  motion  develops  in 
this  loading  geometry.  The  influence  of  lateral  position  within  the  bar  on  the  total  stress  at  a 
range  of  10  in  is  shown  in  Figure  6.23,  both  without  and  with  relative  fluid  motion.  The 
waveforms  are  very  similar,  but  in  the  case  without  motion  slightly  higher  peak  stresses 
develop  toward  the  center  of  the  bar.  In  the  case  with  fluid  motion  these  differences  in  peak 
stress  are  reduced  substantially.  The  influence  of  relative  fluid  motion  on  the  skeleton 
velocities  in  the  uniaxial  stress  loading  is  illustrated  in  Figure  6.23.  As  was  the  case  with 
total  stress,  there  is  very  little  influence  of  the  relative  motion  for  this  loading  condition, 
though  the  relative  fluid  motion  appears  to  damp  out  some  of  the  high  frequency  oscillations 
in  the  calculation  without  relative  motion. 

The  most  dramatic  influence  of  relative  fluid  motion  is  seen  in  the  uniaxial  strain 
calculations.  Figure  6.25  shows  comparisons  of  total  stress  time-histories  at  the  three  ranges 
of  interest  with  and  without  relative  fluid  motion.  The  relative  fluid  motion  dramatically 
alters  the  wave  shapes  by  increasing  the  propagation  velocity,  reducing  the  peak  stress  and 
smearing  both  the  rise  time  and  decay.  As  noted  in  the  previous  discussion,  the  increased 
wavespeed  is  predicted  by  the  closed  form  solution  shown  in  Figure  6.1  (listed  in  Table  6.4). 
Despite  the  rapid  attenuation  of  peak  stress  due  to  the  relative  fluid  motion,  there  appears  to 
be  very  little  reduction  in  total  impulse  transmitted  by  the  stress  wave.  Thus,  the  relative 

fluid  motion  tends  to  spatially  redistribute  the  energy  in  the  wave,  with  very  little  absorption 
of  energy  due  to  the  relative  fluid  motion. 

The  total  stress  waveforms  for  the  calculation  with  relative  fluid  motion  (Figure  6.25) 
are  resolved  into  components  of  effective  axial  stress  and  pore  pressure  in  Figure  6.26.  At  all 
ranges  there  is  a  rapid  rise  in  pore  pressure  conunensurate  with  a  slower  rise  in  effective 
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stress.  This  is  followed  by  a  levelling  off  in  pore  pressure,  then  a  sharp  decline.  The  sharp 
decline  in  pore  pressure  is  accompanied  by  only  a  modest  decline  in  effective  stress  with  the 
total  stress  during  the  later  portion  of  the  waveform  composed  entirely  of  effective  stress. 

The  waveform  comparisons  for  axial  skeleton  velocity  shown  in  Figure  6.27  are  very 
similar  to  those  for  total  stress  shown  in  Figure  6.25.  The  integration  of  the  velocity 
waveforms  in  Figure  6.28  shows  that  the  late  time  displacements  at  all  ranges  are  nearly 
identical  between  the  two  cases  confirming  the  conclusion  that  little  energy  dissipation  occurs 
in  the  case  with  pore  fluid  flow  and  that  the  changes  in  waveform  are  a  result  of  the  relative 
motion  itself  rather  than  due  to  energy  dissipation  mechanisms.  While  the  late  time 
displacements  are  about  equal,  the  rise  lu  peak  displacement  is  much  more  rapid  in  the  case 
of  no  relative  fluid  motion  because  of  the  large  amount  of  dispersion  in  the  calculations  with 
fluid  motion. 

The  motion  of  the  pore  fluid  relative  to  the  porous  skeleton  is  shown  in  Figure  6.29  at 
each  of  the  three  ranges  of  interest.  The  absolute  displacement  of  the  pore  fluid  is  much 
greater  than  that  of  the  skeleton  at  early  time  due  to  the  surge  of  pore  fluid  ahead  of  the 
skeleton.  Eventually,  the  skeleton  overtakes  and  passes  the  pore  fluid,  followed  by  a  very 
gradual  return  of  the  pore  fluid  to  its  original  position  in  the  skeleton  at  very  late  times. 

The  fluid  confined  pressure  bar  calculation  shows  many  of  the  same  characteristics  as 
the  uniaxial  strain  calculation,  but  with  less  relative  fluid  motion  and  dispersion  in  the  case 
where  relative  fluid  motion  is  allowed.  Comparison  of  the  fluid  confined  total  stresses  with 
and  without  relative  motion  is  shown  in  Figure  6.30.  The  case  with  relative  motion  has  a 
faster  wavespeed,  more  peak  attenuation  and  more  dispersion,  all  trends  observed  in  the 
uniaxial  strain  comparisons.  Figure  6.31  demonstrates  that  the  fluid  confined  loading  is  very 
uniform  over  the  entire  width  of  the  bar,  both  with  and  without  relative  motion.  Figure  6.32 
shows  that  the  uniaxial  velocity  time-histories  are  very  similar  to  the  total  stresses.  Finally, 
the  fluid  displacement  and  skeleton  displacement  time  histories  for  the  case  with  relative  fluid 
motion  are  compared  in  Figure  6.33.  The  response  is  again  similar  to  that  for  the  uniaxial 
strain  case  in  that  the  fluid  initially  surges  ahead  of  the  skeleton,  the  skeleton  overtakes  and 


passes  the  fluid,  then  finally  the  two  return  to  their  initial  state  of  relative  equilibrium. 

6.6  CONCLUSIONS  FROM  PARAMETRIC  CALCULATIONS 

The  set  of  calculations  discussed  in  this  section  were  of  value  in  plaiming  and 
designing  the  dynamic  pressure  bar  experiments.  Principal  conclusions  from  the  analysis  of 
these  calculations  include  the  following: 

1.  Pore  fluid  within  an  unconfined  bar  or  within  a  bar  confined  at  a  uniform 
pressure  (such  as  by  confining  the  bar  within  a  pressurized  gas)  has  very  little 
influence  on  the  dynamic  response  of  the  bar.  Because  of  the  immediate  relief 
offered  by  waves  propagating  inward  from  the  sides  of  the  bar,  there  is  almost 
no  buildup  of  pore  fluid  pressure  and  only  very  minor  relative  fluid  motion 
developed.  Neither  configuration  was  considered  satisfactory  for  the 
experimental  portion  of  this  study. 

2.  The  dynamic  uniaxial  strain  loading  condition,  in  which  the  bar  is  perfectly 
restrained  to  zero  expansion  or  contraction  in  the  lateral  direction,  shows 
dramatic  influence  of  relative  pore  fluid  motion  on  both  the  wavespeed  and 
shape  of  the  propagating  wave.  Allowing  pore  fluid  motion,  as  governed  by 
the  stipulated  porosity  and  permeability  of  the  bar,  results  in  increased 
wavespeed,  rapid  attenuation  of  peak  velocities  and  stresses,  and  rapid 
dispersion  in  the  waveform.  These  changes  occur  in  conjunction  with  the 
relative  fluid  motion,  but  are  not  a  result  of  the  minor  amount  of  energy 
dissipation  due  to  fluid  friction.  During  the  dynamic  loading  with  pore  fluid 
motion,  the  pore  fluid  initially  surges  well  ahead  of  the  porous  skeleton, 
followed  by  the  skeleton  overtaking  and  passing  the  pore  fluid  and  finally  the 
fluid  and  skeleton  returning  to  their  initial  relative  positions  very  late  in  the 
loading.  Though  the  uniaxial  strain  loading  is  an  ideal  one  for  studying  the 
influence  of  relative  pore  fluid  motion,  in  practice  it  is  nearly  impossible  to 
satisfactorily  confine  a  pressure  bar  under  such  conditions. 
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3.  Connnement  of  the  pressure  bar  within  a  thin  annulus  of  liquid  confuted  within 
a  steel  vessel  around  the  bar  appears  to  provide  a  satisfactory  nseans  of 
replicating  all  of  the  characteristics  of  the  uniaxial  strain  loading.  Even  though 
the  bar  is  not  perfectly  confined,  confinement  is  sufficient  to  generate  all  of  the 
phenomena  observed  above  in  the  uniaxial  strain  case.  However,  these 
phenomena  are  somewhat  less  prominent  in  the  case  of  the  fluid  confined  bar, 
though  they  should  be  easily  detected  and  measured  in  our  test  apparatus. 


6-19 


Table  6.1.  Summary  of  parametric  calculations. 


Confined  Condition 

Permeable  Porous  Bar 

Impermeable  Porous  Bar 

Uniaxial  Stress  -  Constant 

Confining  Pressure 

WB-U-SR 

WB-U-SN 

Uniaxial  Strain 

WB-C-SR 

WB-C-SN  1 

Liquid  Confinement, 

Small  Annulus 

WB-F-SR 

WB-F-SN  1 
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Table  6.2.  Properties  of  porous  skeleton  and  pore  fluid  used  in  parametric 
calculations. 


Porous  skeleton  (316  L  Stainless  Steel,  SIKA'R-lOOj  | 

Skeleton: 

Porosity: 

n  =  0.465 

Permeability: 

k  =  0.01542  in/sec  (kerosene) 

Drained  Bulk  Modulus: 

K,  =  0.26  X  10*psi 

Poisson’s  Ratio: 

\>  =0.2 

Solid  Grain: 

Specific  Gravity: 

00 

II 

«n 

O 

Bulk  Modulus: 

K,  =  16.67  X  10^  psi 

Pore  Fluid  (Kerosene) 

Specific  Gravity: 

G,  =  0.81 

Bulk  Modulus: 

Kf  =  0.188  X  16*  psi 
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Table  6.4.  Wavespeeds  in  saturated,  porous  bar  determined  from 
parametric  calculations. 


CALCULATION  SET 

RELATIVE 
FLUID  MOTION 
PREVENTED 
(m/s) 

.  I 

RELATIVE 

FLUID  MOTION 
AIXOWED 
(m/s) 

Uniaxial  Stress  (constant  confining  pressure) 

835 

845 

Uniaxial  Strain 

1170 

1375  1 

Liquid  Confinement  (small  annulus) 

1210 

1290  1 
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Damping  (in 


Damping  (in 
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I 

Figure  6.4.  Photograph  of  the  wave  propagation  bar  looking  down  axis.  ^ 

I 

I 
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Figure  6.5.  Photograph  of  the  upstream  end  of  the  incident  bar  showing 
the  axial  load  reaction  frame. 
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Figure  6.6 


Photograph  of  the  incidei 
of  the  pressure  vessel. 
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Photograph  of  the  downstrea.n  end  of  the  pressure  vessel 
showr'v  th<=  electrical  penetrations  and  the  connection 
to  the  specimen  pore  space. 


Stress  (psi) 


Time  (s) 


Figure  6.8.  Stress  waves  measured  in  two  different  locations  on  the 

incident  bar  of  the  wave  propagation  device  due  to  impact 
of  a  20-inch  striker  bar. 
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calculation. 
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10.  Schematic  section  view  of  uniaxial  strain  calculation. 
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Figure  6.11.  Schematic  section  view  of  fluid  confined  calculation. 
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Figure  6.12a.  Influence  of  confinement  condition  on  total  stress  waveforms  with  no  relative  fluid 
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Figure  6.13b.  Influence  of  confinement  condition  on  pore  pressure  with  no  relative 
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Influence  of  confinement  condition  on  pore  pressure  with  no  relative  fluid 
motion  at  Range  C  (15  in.). 
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Figure  6.14c.  Influence  of  confinement  condition  on  particle  velocity  with  no  relative  fluid 
motion  at  Range  C  (15  In.). 
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Figure  6.15a.  Influence  of  confinement  condition  on  particle  displacement  with  no  relative  fluid 
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Figure  6.1bc.  Influence  of  confinement  condition  on  particle  displacement  with  no  relative  fluid 
motion  at  Range  C  (15  in.) 
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Figure  6.16a.  Influence  of  confinement  condition  on  total  stress  waveforms  with  relative  fluid 
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Figure  6.16c.  Influence  of  confinement  condition  on  total  stress  waveforms  with  relative  fluid 
motion  at  Range  C  (15  in.). 
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Figure  6.18a.  Influence  of  confinement  condition  o 
fluid  motion  at  Range  A  (5  in.). 
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Figure  6.18b.  Influence  of  confinement  condition  on  skeleton  particle  veloc 
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Figure  6.18c.  Influence  of  confinement  condition  on  skeleton  particle  velocity  with  relative 
fluid  motion  at  Range  C  (15  In.). 
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Figure  6.19a.  Influence  of  confinement  condition  on  skeleton  particle  displacement  with  relative 
fluid  motion  at  Range  A  (5  in.) 
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Figure  6.19c.  Influence  of  confinement  condition  on  skeleton  particle  displacement  witli  relative 
fluid  motion  at  Range  C  (15  in.). 
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Figure  6.20b.  Influence  of  confinement  condition  on  axial  relative  fluid  velocity  at  Range 
(10  in. ) . 
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Figure  6.20c.  Influence  of  confinement  condition  on  axial  relative  fluid  vejocity  at  Range 
(15  in.). 
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Figure  6.21a.  Influence  of  confinement  condition  on  axial  fluid  displacement  at  Range  A  (5  In 
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Figure  6.22a.  Influence  of  relative  fluid  motion  on  total  stress  lor  uniaxial  stress 
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Figure  6«233«  Influence  of  laternl  position  on  total  stress  for  uniaxial  stress  Joadlng  without 
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Figure  6.24c.  Influence'  of  relative  fluid  motion  on  .skeleton  velocity  for  nni;i 
at  Range  C  (15  in.). 
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Figure  6.25a.  Influence 
loading  a 
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Figure  6.25b.  Influence  of 
loading  at  Ra 
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Figure  6.25c.  Influence  of  relative  fluid  motion  on  total  axial  stress  for  uniaxial  strain 
loading  at  Range  C  (15  in.). 
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Figure  6.26c.  Composite  stress  time-histories  for  uniaxial  loading  with  fluid  motion  at  Range 
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Figure  6.28b.  Influence  of  relative 
loading  at  Range  B  (1 
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Figure  6.28c.  Influence  of  relative  fluid  motion  on  skeleton  displacement  for  uniaxial  strain 
loading  at  Range  C  (15  in.). 
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Figure  6.29a.  Comparison  of  skeleton  and  fluid  displacements  for  uniaxial  strain  loading  at 


Figure  6.29b.  Comparison  of  skeleton  and  fluid  displacements  for  uniaxial  strain  loading  at 
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Figure  6.29c.  Comparison  of  skeleton  and  fluid  displacements  for  uniaxial  strain  loading  at 
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Figure  6.32a.  Influence  jf  relative  fluid  motion  on  skeleton  velocity  for  fluid  confined  loa 
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Figure  6«32c.  Influence  of  relative  fluiil  mo tion  on  skeleton  velticity  lt>r  Iluid  confined  loading 
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Figure  6.33b.  Influence  of  relative  fluid  motion  on  skeleton  displacement  for  fluid  confined 
loading  at  Range  B  (10  in.). 
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Figure  6.33c.  Influence  of  relative  fluid  motion  on  skeleton  displacement  for  fluid  confined 


SECTION  7 


WAVE  PROPAGATION  BAR  TESTS  AND  NUMERICAL  SIMULATIONS 

This  section  presents  the  results  of  a  series  of  tests  that  was  performed  using  the  test 
apparatus  described  in  Section  6  to  investigate  wave  propagation  through  a  fluid-saturated  porous 
medium.  Tests  were  performed  on  a  specimen  of  sintered  stainless  steel  under  both  dry  and 
fluid-saturated  conditions.  This  section  describes  the  porous  stainless  steel  specimen,  including 
mechanical  and  fluid  flow  propenies,  and  presents  the  results  of  the  wave  propagation  tests.  Also 
presented  in  this  section  are  the  results  of  a  numerical  simulation,  performed  with  the  finite 
element  program,  MPDAP,  of  one  wave  propagation  bar  test  on  a  saturated  specimen. 


7.1  WAVE  PROPAGATION  BAR  TEST  DESCRIPTION 

The  Wave  Propagation  Bar  (WPB)  apparatus  used  to  conduct  the  laboratory  experiments 
is  described  in  Section  6.1,  and  depicted  in  Figures  6.2  through  6.7.  The  subsection  describes 
the  sintered  stainless  steel  specimen  and  the  procedures  and  instrumentation  used  to  conduct  the 
tests. 


7.1.1  Porous  Stainless  Steel  Specimen 

All  of  the  WPB  tests  reported  in  this  section  were  performed  on  a  specimen  of  porous 
stainless  steel.  This  material  was  selected  because  it  has  several  advantages  over  a  geologic 
material  for  studies  of  the  fundamental  mechanics  of  wave  propagation.  The  metal  skeleton  is 
more  nearly  linear  elastic  than  any  sand  or  porous  rock.  It  is  very  important  that  the  properties 
of  the  constituent  materials  be  well  characterized  so  that  the  inherent  complexity  of  the  wave 
mechanics  is  not  fur*ner  confused  by  unknown  or  nonlinear  skeleton  response.  Porous  metal  is 
formed  by  processing  granular  material  at  high  temperature  and  pressure  to  form  a  solid  with 
void  spaces  similar  to  a  porous  rock.  Since  this  is  a  controlled  process,  the  properties  of  the 
resulting  material  can  be  specified  within  some  range.  This,  along  with  the  fact  that  the  skeleton 
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is  metal,  makes  it  possible  to  obtain  a  material  that  has  both  relatively  high  permeability  and  high 
skeleton  stiffness,  a  combination  that  was  expected  to  accentuate  the  phenomenology  of  interest 
Blouin,  et.al.  (1990)  listed  fifteen  grades  of  porous  stainless  steel  that  are  commercially  available, 
covering  a  broad  range  of  porosity  and  permeability.  Since  the  elastic  modulus  is  apparently  not 
of  interest  in  the  normal  commercial  uses  of  the  material,  no  data  on  compressibility  properties 
were  available  at  the  time  of  the  initial  design.  The  original  equipment  design  and  calculations 
of  specimen  response  were  based  on  estimates  of  the  mechanical  properties  of  the  various  grades 
of  porous  stainless  steel,  as  described  in  Section  6. 

Two  bars  of  porous  stainless  steel  were  fabricated  for  use  as  test  specimens  by  Newmet 
Krebsoge  of  Terry  ville,  CT.  The  nominal  dimensions  of  the  bars  were  2  inches  diameter  and  24 
inches  long.  One  bar  was  cut  into  shorter  lengths  for  mechanical  property  and  fluid  flow  testing, 
and  the  other  was  used  at  almost  full  length  for  WPB  tests.  In  order  to  avoid  closing  the  pore 
spaces  on  the  cut  as  would  happen  with  conventional  cutting  techniques,  the  bars  were  cut  using 
a  modem  machining  technique  known  as  wire  electric  discharge  machining. 

The  bars  were  specified  to  be  Sika  R  100.  However,  since  the  material  is  not  normally 
manufactured  in  such  large  sizes,  some  variation  from  the  normal  R  100  specifications  was 
expected.  The  nominal  porosity  of  the  Sika  R  100  based  on  the  specified  density  is  0.465,  and 
it  was  estimated  to  have  an  elastic  modulus  of  4.7x10’  psi  and  Poisson’s  ratio  of  0.2.  In  order 
to  determine  the  actual  properties  of  the  test  specimens,  an  unconflned  compression  test  was 
performed  on  a  cylindrical  specimen  of  the  porous  stainless  steel  (approx.  2  in.  diameter  by  4  in. 
long).  Plots  of  the  resulting  axial  and  radial  strain  against  axial  stress  are  presented  in  Figure 
7.1.  From  the  data  in  Figure  7.1  the  elastic  modulus  and  Poisson’s  ratio  were  found  to  be  3.6 
X  10®  psi  and  0.25,  respectively.  Thus,  the  stiffness  of  the  porous  bars  as  fabricated  was  nearly 
and  order  of  magnitude  greater  than  originally  estimated.  Confirmation  of  the  unconfmed 
compression  test  result  is  provided  by  data  acquired  while  preparing  the  full  length  specimen  for 
wave  propagation  testing.  The  strain  gage  circuits  on  the  porous  specimen  and  axial  load 
measurement  gages  on  the  incident  bar  were  monitored  during  initial  static  axial  preloading  of 
the  specimen.  The  elastic  modulus  computed  from  those  measurements  was  found  to  be 
essentially  identical  to  the  value  determined  from  the  unconfmed  compression  test  on  the 
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specimen  cut  from  the  other  bar  of  porous  material.  The  physical  and  mechanical  properties  of 
the  porous  specimen  are  summarized  in  Table  7.1. 


The  results  of  a  series  of  permeability  tests  on  a  specimen  of  the  porous  stainless  steel 
are  presented  in  Figure  7.2.  A  straight  line  was  fit  using  the  complete  data  set,  as  depicted  in 
the  figure.  The  intercept  and  slope  of  the  line  were  used  to  define  flow  coefficients  a  and  b, 
respectively,  as  defined  in  Section  4.  The  permeability  is  lower  by  approximately  a  factor  of  five 
than  the  value  given  by  the  manufacturer,  reported  by  Blouin,  et  al.,  (1990),  and  used  in  the 
parametric  calculations  reported  in  Section  6.2. 

The  expected  wavespeed  and  damping  for  compressional  waves  of  the  first  and  second 
kinds  in  kerosene  saturated  porous  stainless  steel  with  the  measured  material  properties  were 
computed  as  a  function  of  the  frequency  of  the  excitation  pulse  with  the  computer  program 
TWAVE,  as  described  in  Section  6.  The  results  of  those  calculations  are  presented  graphically 
in  Figure  7.3.  Table  7.2  presents  the  wavespeeds  of  the  saturated  material  and  its  various 
constituents,  computed  from  the  properties  listed  in  Table  7.1,  using  the  procedures  described  in 
Section  6.6.  Comparison  of  the  values  in  Table  7.2  with  the  corresponding  quantities  in  Table 
6.3  reveals  that  the  wavespeeds  computed  based  on  the  initial  estimates  of  specimen  stiffness  are 
significantly  in  error. 

The  transition  in  the  wave  of  the  first  kind  occurs  in  the  range  of  approximately  10,000 
to  100,000  rad/s.  Thus,  the  range  of  frequencies  possible  with  the  wave  propagation  apparatus 
as  fabricated  (15,700  to  62,8(X)  rad/s)  falls  within  the  transition  region  for  the  specimen  material. 


7.1.2  Test  Procedures 

The  specimen  was  21.3  inches  long,  and  its  cylindrical  surface  was  rough  due  to  the 
fabrication  process,  having  a  diameter  of  1.90  ±  0.02  inches.  Prior  to  installation  of  the  specimen 
in  the  test  apparatus,  two  foil  strain  gages,  one  each  in  the  axial  and  circumferential  directions, 
were  bonded  to  the  bar  at  specially  prepared  locations  on  its  surface  at  distances  of  5,  10,  and 
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IS  inches  from  the  upstream  end  of  the  specimen  and  a  pore  pressure  gage  was  also  placed  on 
the  surface  of  the  specimen  at  the  5-inch  station,  as  described  in  Section  7.1.3.  The  specimen 
was  then  placed  inside  a  heat-shrinkable  polyolefin  tube  and  installed  between  the  incident  bar 
and  transmitter  bar  as  illustrated  in  Figure  6.2.  The  tubing  was  shrunk  to  a  tight  fit  around  the 
specimen.  The  tubing  also  extended  approximately  1  inch  onto  the  bars  on  either  end  and  was 
sealed  to  the  incident  and  transmitter  bars  with  epoxy  adhesive  and  wire  clamps.  At  this  point, 
the  pressure  vessel  was  closed,  and  all  of  the  instrumentation  was  connected  and  checked.  Wave 
propagation  measurements  were  made  on  the  dry  specimen  with  no  confining  pressure.  During 
these  tests,  approximately  350  psi  average  axial  stress  was  applied  to  the  specimen  to  tightly  seat 
the  various  components  of  the  WPB  assembly. 

Upon  completion  of  the  dry  tests,  the  pressure  vessel  was  evacuated  and  then  filled  with 
de-aired  kerosene.  Access  to  the  pore  space  of  the  specimen  was  by  means  of  an  axial  hole 
drilled  in  the  transmitter  bar  and  an  intersecting  radial  hole  outside  the  pressure  vessel.  In 
preparation  for  saturation,  the  pore  space  was  evacuated.  It  was  then  filled  with  de-aired 
kerosene.  A  quantity  of  kerosene  slightly  larger  than  the  void  space  of  the  specimen  was 
introduced,  resulting  in  a  layer  of  fluid  between  the  specimen  and  the  polyolefin  membrane.  The 
confining  fluid  in  the  pressure  vessel  surrounding  the  specimen  was  then  pressurized  to  1000  psi. 
Axial  load  equivalent  to  1000  psi  average  axial  stress  was  also  applied  to  the  specimen.  Because 
of  the  excess  fluid  surrounding  the  specimen,  this  also  had  the  effect  of  pressurizing  the  pore 
fluid.  This  tends  to  force  any  remaining  pore  air  into  solution.  The  pore  fluid  valve  was  then 
opened  momentarily,  allowing  the  excess  pore  fluid  to  drain.  The  quantity  of  fluid  drained  was 
manually  controlled  to  obtain  the  desired  initial  pore  pressure  of  500  psi.  These  initial  stress 
conditions  of  1000  psi  confining  pressure  and  axial  stress,  and 

500  psi  pore  pressure  were  maintained  throughout  the  wave  propagation  tests  on  saturated 
specimens. 
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7.U  Instrumentation 


During  each  test,  a  high-speed  digital  data  acquisition  system  was  used  to  record  eight 
channels  of  instrumentation,  consisting  of  strain  gages  on  the  porous  stainless  steel  specimen, 
pore  pressure  at  the  surface  of  the  specimen,  and  a  strain  gage  bridge  on  the  incident  bar  to 
measure  the  loading  applied  to  the  test  specimen.  The  locations  of  the  instruments  on  the  porous 
specimen  and  on  the  incident  bar  are  shown  in  Figure  7.4. 

Standard  foil  strain  gages  were  bonded  to  the  stainless  steel  were  used  to  monitor  skeleton 
deformation.  As  manufactured,  the  surface  of  the  porous  stainless  steel  is  much  too  rough  to  be 
strain  gaged.  However,  experimentation  in  the  laboratory  showed  that  by  grinding  the  surface, 
it  is  readily  possible  to  obtain  a  surface  that  can  be  strain  gaged.  The  grinding  not  only  smooths 
the  roughness  but  also  smears  some  material  into  the  void  spaces,  creating  a  surface  with  an 
appropriate  texture  for  strain  gaging.  Axial  and  radial  strains  were  measured  at  there  stations. 
The  three  stations,  designated  A,  B,  and  C  were  located  at  distances  of  5,  10,  and  15  inches, 
respectively,  from  the  incident  end  of  the  test  specimen. 

One  pore  pressure  measurement  was  made  at  the  surface  of  the  porous  bar  at  Station  A, 
located  5  inches  from  the  upstream  end  of  the  specimen.  A  special  instrumentation  package  was 
devised  for  that  purpose.  An  ENTRAN  Model  EPL6-125-5000X  pressure  transducer  (0. 125  inch 
diameter  and  0.07  inch  thick)  was  installed  under  a  specially  designed  aluminum  shroud  which 
distributed  the  confining  pressure  load  around  the  pressure  transducer  so  that  it  was  exposed  only 
to  pore  pressure. 

The  load  applied  to  the  test  specimen  was  measured  with  a  full  strain  .gage  bridge  on  the 
incident  bar  external  to  the  pressure  vessel.  The  foil  strain  gages  were  bonded  to  the 
incident  bar  37.5  inches  upstream  from  the  contact  with  the  test  specimen.  They  were  arranged 
and  wired  in  the  conventional  manner  for  a  colunrn  load  cell. 

All  channels  were  digitized  at  a  sampling  interval  of  1  p.s.  In  each  test,  the  start  time  for 
data  acquisition  was  initiated  when  the  striker  bar  made  contact  with  a  switch  in  advance  of 
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impacting  the  incident  bar.  Due  to  variations  in  switch  positioning  and  bar  velocity,  this  is  not 
a  reliable  zero-time  indication.  In  all  of  the  test  data  plots,  zero  time  has  been  set  to 
approximately  the  time  the  stress  wave  arrived  at  the  interface  between  the  incident  bar  and  the 
test  specimen.  This  was  done  by  choosing  a  point  55  ps  before  arrival  of  the  wave  at  strain  gage 
Staion  A.  It  is  estimated  that  the  difference  between  the  zero  time  selected  in  this  manner  and 
the  actual  arrival  of  the  stress  wave  is  less  than  5  ps.  All  of  the  traces  for  each  test  are  plotted 
with  respect  to  the  same  time  reference. 


7.2  WAVE  PROPAGATION  BAR  TEST  RESULTS 

This  subsection  presents  the  results  of  the  WPB  tests  on  a  porous  stainless  steel  specimen 
described  in  Section  7.1.  Table  7.3  presents  a  summary  of  all  the  tests  performed.  The  results 
of  tests  on  the  porous  specimen  m  the  dry  condition  are  presented  in  Section  7.2.1,  and  Section 

7.2.2  presents  the  tests  performed  with  the  specimen  saturated. 

7.2.1  Wave  Propagation  Tests  on  Dry  Specimen 

As  indicated  in  Table  7.3,  tests  on  dry  porous  stainless  steel  specimens  were  performed 
with  three  different  length  striker  bars,  5-,  10-,  and  20-inch.  Three  nominally  identical  tests  were 
performed  with  each  striker  bar.  There  were  no  significant  differences  between  tests  with  a  given 
bar  type,  and  a  typical  test  of  each  type  was  selected  for  further  analysis.  Figures  7.5  through 
7.7  present  the  results  of  three  typical  tests  on  dry  specimens,  one  with  each  striker  bar  length. 
Each  of  Figures  7.5  through  7.7  has  three  parts,  designated  a,  b,  and  c.  Part  a  of  each  figure 
presents  the  input  force  time  history  measured  on  incident  bar  for  that  test.  Parts  b  and  c  of 
Figures  7.5  through  7.7  present  axial  and  radial  strains,  respectively,  on  the  test  specimen.  Each 
of  the  strain  plots  shows  three  strain  histories,  corresponding  to  the  three  strain  measurement 
stations  as  defined  above  and  illustrated  in  Figure  7.4. 

By  examining  parts  a  of  Figures  7.5  through  7.7,  it  can  be  seen  that,  as  expected, 
durations  of  the  applied  stress  pulses  are  proportional  to  their  respective  striker  bar  lengths.  The 
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S-,  10-,  and  20-iach  striker  bars  produced  pulse  durations,  as  measured  from  the  beginning  of  the 
rise  to  the  onset  of  unloading,  are  50,  100,  and  200  ps,  respectively.  All  three  tests  exhibit 
wavespeeds  of  about  90,000  in./s  (2300  m/s)  when  determined  from  first  arrivals  and 
approximately  80,000  in./s  (2000  m/s)  when  measured  from  the  mid-points  of  the  initial  rise 
portions  of  the  stress  wave  signals.  The  wavespeed  determined  from  first  arrivals  is  in  good 
agreement  with  the  rod  (uniaxial  stress)  wavespeed  listed  in  Table  7.2.  This  result  differs  from 
the  numerical  simulation  results  presented  in  Section  6,  in  which  the  wavespeed  at  mid-height 
of  the  first  rise  to  peak  was  in  reasonable  agreement  with  the  computed  uniaxial  strain 
wavespeed.  Recall  that  the  actual  elastic  modulus  of  the  porous  stainless  steel  is  almost  a  factor 
of  eight  larger  than  the  elastic  modulus  used  in  the  calculations  presented  in  Section  6.  The 
difference  in  wavespeeds  is  an  indication  that  the  fluid  confinement  is  not  as  effective  in 
maintaining  a  uniaxial  strain  condition  with  the  real  material  as  it  was  computed  to  be  with  the 
erroneously  soft  porous  skeleton. 

The  reflections  from  the  downstream  end  of  the  specimen  are  evident  in  the  strain  data 
presented  in  Figures  ”^.5  through  7.7  (parts  b  and  c).  The  reflection  arrives  first  at  Station  C,  at 
approximately  350  ps.  It  is  followed  by  arrivals  at  Stations  B  and  A  at  times  of  approximately 
400  and  450  ps,  respectively. 

7.2.2  Wave  Propagation  Tests  on  Saturated  Specimen 

Wave  propagation  bar  tests  were  performed  on  the  same  porous  stainless  steel  specimen 
fully  saturated  with  kerosene,  as  described  in  Section  7.1.  Those  tests  are  summarized  in  Table 
7.3.  In  addition  to  tests  with  the  5-,  10-,  and  20-inch  striker  bars,  a  test  was  performed  on  the 
saturated  specimen  with  a  1/8-inch  thick  rubber  pad  attached  to  the  front  of  the  20-inch  striker 
bar  to  produce  an  input  stress  pulse  with  a  longer  rise  time  to  peak.  Figures  7.8  through  7.11 
present  data  plots  for  one  representative  test  with  each  type  of  excitation.  Figures  7.8  through 
7.11  correspond  respectively  to  excitation  by  5-inch  striker,  10-inch  striker,  20-inch  striker,  and 
20-inch  striker  with  rubber  pad.  As  with  the  dry  specimens,  each  figure  has  parts  a,  b,  and  c 
showing  the  input  force  pulse,  axial  strains,  and  circumferential  strains,  respectively.  In  addition, 
each  figure  for  the  saturated  tests  has  a  part  d  which  presents  the  pore  pressure  time  history  at 
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Station  A,  5  inches  from  the  interface  with  the  incident  bar. 


The  wavespeed  in  the  saturated  porous  bar,  as  determined  from  the  first  arrivals  is  90,000 
in./s  (2300  m/s).  When  determined  by  the  mid-points  of  the  initial  rise  portions  of  the  waves  at 
the  first  and  last  strain  gage  stations,  the  wavespeed  is  80,000  in./s  (2000  m/s).  These  values  are 
the  same  as  the  corresponding  measurements  for  the  dry  porous  bar.  This  is  consistent  with  the 
calculated  wavespeeds  presented  in  Table  7.2,  which  show  no  significant  difference  in  computed 
wavespeeds  between  the  dry  and  saturated  cases.  This  lack,  of  sensitivity  of  wavespeed  to  the 
degree  of  saturation  results  from  the  dominant  stiffness  of  the  material  skeleton. 

One  measure  of  the  consistency  between  the  experimental  results  and  the  analytical 
understanding  of  the  physics  of  the  problem  is  the  amount  of  the  load  carried  by  the  water 
relative  to  the  porous  skeleton.  Blouin  and  Kim  (1984)  present  the  following  expression  for  the 
ratio,  P„,,  of  effective  stress  to  pore  pressure  under  uniaxial  strain  loading  given  by: 

Kl  .  g,  g)  -  M,  K,  K,  -  K,  K,  K, 

jf.  a:,  («,  -  at; 


where:  =  mixture  modulus  =  (K^  K^)  /[K^  +  n(Kg  -  K^)] 

Kg,  Kn,,  K,  =  bulk  modulus  of  grains,  fluid,  skeleton 

M,  =  constrained  modulus  of  skeleton  =  3K,(l-v)/(l-hv) 

V  =  Poisson’s  ratio  of  skeleton 
n  =  porosity  of  skeleton 

Using  the  properties  for  the  porous  stainless  steel  skeleton,  solid  grains,  and  kerosene  given  in 
Table  7.1,  the  resulting  value  of  is  9.75.  For  comparison,  consider  the  test  data  from  a 
saturated  specimen  impacted  by  a  10-inch  striker  bar  which  are  shown  in  Figures  7.9b-7.9d.  The 
peak  in  the  pore  pressure  history  occurs  at  approximately  120  p,s.  The  initial  peaks  in  both 
skeleton  strains  at  Station  A  occur  at  the  same  time.  Using  the  axial  and  circumferential  strains 
at  120  jis  and  the  elastic  constants  given  in  Table  7.1,  the  axial  stress  in  the  stainless  steel 
skeleton  is  found  to  be  approximately  3500  psi.  At  the  same  time  the  pore  pressure  was 
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approximately  345  psi.  The  ratio  of  effective  stress  in  the  skeleton  as  derived  from  strain 
measurements  to  the  measured  dynamic  pore  pressure  at  the  same  time  is  10.1.  This  quantity 
is  comparable,  although  not  exactly  the  same,  as  the  derived  from  theory.  The  ratio,  P„,  is 
derived  assuming  a  uniaxial  strain  condition.  While  this  is  known  not  to  be  the  case  in  the 
experiment,  the  close  match  between  measured  and  calculated  values  is  a  clear  indication  that 
the  porous  skeleton  was  fully  saturated  and  that  the  partitioning  of  load  between  skeleton  and 
pore  fluid  was  approximately  as  predicted  by  theory. 


7.2.3  Comparison  between  Dry  and  Saturated  Test  Results 

Figures  7.12  through  7.14  present  comparisons  of  strain  records  from  dry  and  saturated 
tests  for  striker  bar  lengths  of  5,  10,  and  20  inches,  respectively.  Parts  a  of  Figures  7.12  through 
7.14  show  axial  strains  and  parts  b  present  circumferential  strains.  The  traces  have  been  shifted, 
as  necessary  so  that  the  initial  arrivals  occur  at  the  same  time.  The  records  from  the  dry  tests 
in  Figures  7.13  and  7.14  were  also  scaled  in  amplitude  by  the  ratios  of  the  loading  amplitudes 
so  that  they  are  approximately  the  same  size  on  the  comparison  plots.  The  loading  amplitudes 
for  the  dry  and  saturated  tests  were  within  2%  of  each  other  for  the  records  shown  in  Figure 
7.12,  and  no  scaling  was  done.  In  all  three  tests,  the  arrivals  of  the  waves  at  Stations  B  and  C 
occurred  at  essentially  the  same  times  in  both  dry  and  saturated  specimens.  This  shows,  as  noted 
in  the  preceding  subsections,  that  the  wavespeeds  in  the  saturated  tests  were  essentially  identical 
to  those  in  the  dry  tests. 

In  the  tests  with  the  5-inch  striker  bar,  shown  in  Figure  7.12,  the  wave  propagating 
through  the  dry  specimen  exhibited  more  dispersion  than  saturated  specimen,  as  evidenced  by  the 
longer  durations  of  the  strain  pulses.  This  trend  does  not  appear  to  carry  over  to  the  tests  with 
the  longer  striker  bars  shown  in  Figures  7.13  and  7.14. 
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13  PrtJMERICAL  SIMULATION  OF  LABORATORY  EXPERIMENT 


To  improve  our  understanding  of  the  laboratory  results  and  to  validate  our  numerical 
models,  a  detailed  finite  element  calculation  was  performed  to  simulate  wave  propagation  test 
number  D31K1.  The  calculation  was  completed  using  the  MPDAP  code  and  was  similar  to  the 
parametric  calculations  discussed  in  Section  6.  The  finite  element  mesh  geometry  and  boundary 
conditions  were  the  same  as  used  in  the  parametric  calculations  of  the  fluid  confined  bar  as 
depicted  in  Figure  6.11.  For  this  simulation,  however,  the  material  properties,  initial  stress 
conditions,  and  input  loading  wave  were  carefully  selected  to  model  the  actual  test  conditions  as 
closely  as  possible. 

The  porous,  stainless  steel  test  specimen  and  the  kerosene  fluid  were  modeled  as  linear 
elastic  materials.  This  is  an  appropriate  assumption  for  the  stress  wave  amplitudes  induced  in 
the  tests.  Properties  for  the  test  specimen  and  pore  fluid  are  summarized  in  Table  7.1.  Elastic 
properties  for  the  porous  steel  skeleton  were  obtained  from  the  unconfirmed  compression  test  on 
a  sample  of  the  specimen  material,  as  described  in  Section  7.1.1  and  shown  in  Figure  7.1.  The 
specific  gravities  and  bulk  moduli  for  the  solid  steel  grains  and  the  kerosene  were  obtained  from 
standard  references  on  material  properties.  Note  that  kerosene  is  used  for  both  the  pore  fluid 
inside  the  porous  skeleton  and  the  confining  fluid.  The  permeability  parameters  for  the  porous- 
steel  test  specimen  are  also  given  in  Table  7.1.  They  were  derived  from  the  test  data  presented 
in  Figure  7.2.  Since  the  porosity  of  the  bar  is  practically  unchanged  during  passage  of  the  stress 
wave,  the  permeability  parameters  were  considered  as  constant. 

The  parametric  calculations  described  in  Section  6.4  were  started  with  a  zero  stress 
condition  in  the  specimen.  In  test  number  D31K1,  an  axial  stress  of  1000  psi  was  applied  to  the 
test  specimen  to  ensure  full  contact  with  the  incident  and  transmitter  bars.  Furthermore,  a 
500  psi  initial  pore  pressure  was  applied  to  the  specimen  to  maintain  saturation  necessary  for 
accurate  pore  pressure  measurements.  Both  the  initial  axial  stress  of  1000  psi  and  pore  pressure 
of  500  psi  were  applied  at  the  start  of  the  numerical  simulation  to  maintain  a  complete  simulation 
of  the  test.  However,  since  only  linear  elastic  materia]  models  were  used,  modeling  the  initial 
stress  condition  should  have  no  effect  on  the  simulation  results. 


To  permit  a  direct  comparison  between  the  simulation  and  test  results,  it  was  necessary 
to  match,  as  closely  as  possible,  the  loading  conditions  in  the  test.  First,  the  loading  was 
modeled  by  specifying  uniform  velocity  time  history  across  the  end  of  die  bar.  The  velocity 
boundary  ensures  that  the  end  of  the  specimen  will  displace  uniformly  in  the  simulation.  This 
is  a  more  accurate  model  of  the  actual  loading  condition  than  the  pressure  loading  used  in  the 
parametric  calculations.  Further,  the  input  loading  was  applied  only  to  the  end  of  the  specimen 
and  not  to  the  confining  fluid  as  was  done  in  the  parametric  calculations  reported  in  Section  6.4. 

To  develop  an  appropriate  boundary  velocity-time  history  to  simulate  the  stress  wave 
loading,  the  input  stress  pulse  was  recorded  using  a  load  cell  on  the  incident  bar.  From  the 
recorded  stress  wave,  seven  time  points  were  used  directly  in  the  load  history  for  the  simulation. 
To  determine  the  appropriate  velocity  magnitude  at  each  point,  a  trial  and  error  procedure  was 
used  to  determine  a  correlation  factor  between  the  recorded  axial  load  and  the  corresponding 
particle  velocity  to  be  used  in  the  simulation.  In  this  manner,  a  seven-point,  velocity  time  history 
was  developed  that  accurately  simulated  the  input  pulse  in  test  number  D31K1.  This  loading 
condition  is  plotted  in  Figure  7.15. 

The  calculated  axial  strains  in  the  test  specimen  are  overlaid  with  the  actual  test  data  from 
D31K1  in  Figure  7.16a,  b,  and  c  for  ranges  of  5  inch,  10  inch,  and  15  inch,  respectively. 
Consideiing  the  complexity  of  the  problem,  the  overall  agreement  between  the  calculation  and 
the  data  is  excellent.  The  wavespeeds  in  the  simulation  are  slightly  higher  than  those  measured 
in  the  experiments.  The  wavespeed  calculated  based  on  initial  arrival  is  95,000  in./s  (2410  m/s), 
as  compared  with  90,000  in./s  (2290  m/s)  measured  in  the  experiments.  Similarly,  at  the  mid¬ 
points  of  the  initial  rise  to  peak  in  the  strain  signals,  the  calculation  shows  a  speed  of  92,000  in./s 
(2340  m/s)  while  the  rise  to  peak  in  the  tests  was  slower,  exhibiting  a  wavespeed  at  the  mid-point 
of  approximately  80,000  in./s  (2030  m/s).  The  amplitudes  of  both  the  initial  peak  and  the  overall 
peak  are  in  quite  good  agreement  between  calculation  and  experiment. 

Figures  7.17a,  b,  and  c  present  comparisons  of  circumferential  strain  records  at  Stations 
A,  B,  and  C,  respectively.  The  numerical  simulation  under-predicted  the  circumferential  strain 
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at  every  station  relative  to  the  test  measurements.  This  is  an  indication  that  the  confinement 
system  was  somewhat  more  effective  than  the  calculatior-  would  suggest. 

In  Figures  7.16a  and  7.16b  is  evidence  of  what  may  be  the  wave  of  the  second  kind.  The 
simulation  clearly  shows  an  initial  peak  in  axial  strain  followed  by  a  second  significant  rise.  At 
Station  A  (Figure  7.16a),  the  second  rise  begins  at  lOS  ps,  which  is  10  |xs  after  the  initial  peak. 
At  Station  B  (Figure  7.16b),  the  second  rise  begins  at  205  |is,  which  is  about  55  (is  after  the 
initial  peak.  This  secondary  strain  pulse  is  propagating  at  a  velocity  of  50,000  in./s  (4170  ft/s 
=  1270  m/s),  which  agrees  very  well  with  the  computed  velocity  of  the  wave  of  the  second  kind 
(4100  ft/s)  at  the  high  frequency  end  of  Figure  7.3b.  At  the  time  where  this  wave  pulse  would 
be  expected  at  Station  C  (approx.  350  |is),  the  reflection  from  the  far  end  of  the  bar  dominates 
the  waveform,  making  it  impossible  to  detect  the  slow  moving  forward  propagating  pulse.  The 
same  phenomenon  is  present  in  the  lab  test  data,  although  somewhat  less  pronounced.  The 
arrival  of  the  secondary  wave  is  easier  to  observe  in  the  circumferential  strain  data  shown  in 
Figures  7.17a  and  7.17b  than  in  the  axial  strain  records.  In  the  test  data,  the  second  wave  is 
propagating  at  about  40,6(X)  in/s  (3380  ft/s  =  1030  m/s).  As  with  the  wave  of  the  first  kind,  this 
is  somewhat  lower  than  the  speed  predicted  by  analysis  for  the  case  where  radial  strains  are  fully 
constrained. 

If  the  secondary  wave  discussed  in  the  previous  paragraph  is  in  fact  the  wave  of  the 
second  kind,  which  is  associated  with  a  fluid  saturated  porous  material,  then  it  should  not  appear 
in  the  tests  performed  on  the  dry  bar.  Attention  is  referred  back  to  Figures  7.14a  and  7.14b.  At 
Station  A,  the  saturated  waveforms  for  both  axial  and  radial  strains  clearly  exhibit  a  secondary 
pulse  that  is  different  than  those  for  the  dry  material.  At  Station  B,  which  is  5  inches  further 
downstream,  the  secondary  wave  is  significantly  attenuated,  but  still  appears  in  the  records  from 
the  saturated  specimen  and  not  those  for  the  dry  condition. 

Figure  7.18  compares  the  measured  and  calculated  dynamic  pore  pressures.  (In  this 
figure,  zero  dynamic  pore  pressure  corresponds  to  the  initial  pore  pressure  of  500  psi).  In  this 
case,  the  calculated  pore  pressures  are  about  one  half  of  the  measured  values.  Only  one  pore 
pressure  measurement  location  was  available  in  this  test,  so  only  the  comparison  in  Figure  7.18 
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is  possible.  The  computed  pore  pressure  responses  at  5,  10,  and  15  inch  ranges  are  compared 
in  Figure  7.19.  At  the  times  of  the  secondary  strain  pulses  discussed  in  the  previous  two 
paragraphs,  simulation  results  exhibit  a  sharp  drop  in  pore  pressure  at  both  locations.  In  the 
measured  strain  records  shown  in  Figure  7.14,  the  arrival  of  the  secondary  wave  occurs  at 
approximately  130  ps.  Consistent  with  the  analysis  results,  the  one  measured  pore  pressure 
record  also  has  a  major  negative-going  excursion  at  that  time.  However,  in  the  simulations,  the 
first  negative-going  pulse  in  the  pore  pressure  corresponds  with  the  arrival  of  the  secondary  wave 
in  the  strain  records.  In  contrast,  in  the  laboratory  test  records,  the  pore  pressure  contains  a  prior 
negative-going  pulse. 

Figures  7.20a  through  7.20c  present  time  histories  of  particle  velocity  of  the  porous 
skeleton  and  of  velocity  of  the  pore  fluid  relative  to  the  skeleton  at  Stations  A,  B,  and  C, 
respectively.  Particularly  at  Stations  A  and  B,  the  skeleton  velocities  exhibit  secondary  pulses 
of  the  same  form  as  the  strain  records.  In  the  apparent  relative  fluid  velocity  records,  both  the 
arrival  of  the  initial  wave  and  the  secondary  wave,  there  are  short  duration  (less  than  50  ps) 
negative  pulses  of  relative  fluid  velocity,  i.e.  the  fluid  momentarily  lags  behind  the  skeleton. 
Unfortunately,  there  are  no  analogous  laboratory  measurements  for  comparison. 

Figures  7.21  and  7.22  compare  the  computed  axial  strains  and  dynamic  pore  pressures  at 
various  locations  within  the  bar  at  a  range  of  10  inches.  Since  measurements  can  only  be 
acquired  at  the  surface  of  the  specimen  in  the  test,  the  finite  element  simulation  can  be  used  to 
study  variations  in  the  stress  wave  with  distance  from  the  surface  of  the  bar.  From  Figures  7.21 
and  7.22,  it  appears  that  a  very  uniform  stress  condition  is  propagating  through  the  bar  in  this 
simulation. 


7.4  DISCUSSION  OF  WAVE  PROPAGATION  BAR  TESTS 

The  research  effort  described  in  this  section  involves  both  laboratory  experiments  and 
numerical  simulations  of  wave  propagation  through  saturated  porous  media.  A  laboratory  test 
apparatus  was  designed  and  assembled  and  a  series  of  tests  was  performed  on  a  specimen  of 
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porous  (sintered)  stainless  steel.  The  tests  were  conducted  with  the  specimen  both  dry  and  fully 
saturated  with  kerosene.  Instrumentation  on  the  specimen  was  used  to  measure  axial  and 
circumferential  strains  at  three  stations  and  pore  fluid  pressure  at  one  station.  Mechanical  waves 
were  induced  in  the  test  specimen  by  impacting  a  hardened  steel  bar  that  was  in  contact  with  the 
specimen  by  second  hardened  steel  bar  called  the  striker  bar.  The  duration  of  the  stress  pulse 
was  determined  by  the  length  of  the  striker  bar,  and  strikers  of  three  different  lengths  were  used 
in  the  tests  to  produce  stress  pulses  of  nominally  50.  100,  and  200  ps  duration. 

The  numerical  simulations  were  performed  with  MPDAP,  a  special  purpose  finite  element 
computer  program  specifically  designed  for  the  analysis  of  fully  or  partially  saturated  porous 
materials.  A  simulation  was  performed  with  the  test  apparatus  and  input  boundary  conditions 
carefully  defined  to  match,  as  closely  as  possible,  one  of  the  laboratory  experiments. 
Comparisons  were  made  between  the  measured  and  computed  response,  and  between  measured 
response  of  the  porous  stainless  steel  specimen  in  the  dry  and  fully  saturated  conditions. 

In  a  saturated  porous  material,  the  stress  resulting  from  an  imposed  load  is  partitioned 
between  the  porous  skeleton  and  the  pore  fluid.  In  the  test  apparatus,  the  stresses  carried  in  by 
the  solid  and  fluid  phases  were  measured  and  the  distribution  of  stresses  between  the  two  was 
found  to  be  consistent  with  theoretical  predictions.  Unfortunately,  because  of  the  very  high 
stiffness  of  the  porous  stainless  steel  specimen  with  respect  to  the  pore  fluid,  approximately  90% 
of  the  load  was  carried  by  the  steel  skeleton.  Thus,  while  the  saturated  specimen  did  behave  as 
expected,  its  response  was,  in  many  ways,  not  very  different  from  that  of  a  dry  specimen. 
However,  using  the  results  of  the  numerical  simulations  for  guidance,  we  were  able  to  find 
evidence  in  the  both  the  strain  and  pore  pressure  records  from  the  laboratory  experiments  of  the 
wave  of  the  second  kind  that  has  been  theoretically  predicted  to  exist  in  saturated  porous 
materials. 

We  believe  that  our  ability  to  resolve  the  phenomena  of  interest  in  the  test  data  was 
limited  because  of  the  very  high  stiffness  of  the  porous  skeleton  that  we  inadvertently  chose  for 
the  experiment.  In  future  work  it  is  recommended  that  softer  (less  stiff)  skeleton  materials  be 
used  in  order  to  accentuate  the  effect  of  the  pore  water  on  the  overall  response  of  the  saturated 
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porous  material.  Further,  it  would  be  very  desirable  to  devise  an  instrumentation  schenre  that 
would  measure  the  velocity  of  the  pore  fluid  relative  to  the  porous  skeleton. 


745 


Table  7.1.  Properties  of  test  specimen  used  in  numerical 
simulation  of  laboratory  test 


Porous  Stainless  Steel 

Model:  linear  elastic 

porosity 

n 

0.42 

Poisson’s  ratio 

V 

0.25 

Young’s  modulus 

E 

3.6  X  10*  psi  II 

bulk  modulus 

K 

2.4  X  10*  psi 

shear  modulus 

G 

1.44  X  10*  psi 

constrained  modulus 

M 

4.32  X  10*  psi 

Solid  Stainless  Steel  Grains 

Model:  linear  elastic 

specific  gravity 

G, 

7.8 

bulk  modulus 

_ h _ 

16.67  X  10*  psi 

n 

Kerosene  Pore  Fluid  and  Confining  Fluid 

Model:  linear  elastic 

specific  gravity 

G, 

0.81  1 

bulk  modulus 

K, 

0.188  X  10*  psi  1 

viscosity 

It 

2.65  X  10  ’  Ib-s/in’  | 

Permeability  of  Skeleton 

Permeability  coefficient 

a 

7.8  Ib-s/in^ 

Permeability  coefficient 

b 

2.3  Ib-sW 

Permeability  coefficient 

P  =  Pf/b 

3.4  X  10-*  in^ 

Permeability  coefficient 

a  =  p/a 

4.1  X  10  *  in 

7-16 


MATERIAL  CONDITION  MODULUS  SPECIFIC  WAVESPEED 

(psi)  GRAMTY  (m/s) 


Table  7.3.  Wave  propagation  bar  test  summary 
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Apparent  Flow  Velocity  (in/s) 

igure  7.2.  Permeability  teat  results  tor  a  spui;  i  (iii*u  ol  t  lu-  ixirdu.-,  st  <i  i  n  I  I'.ss  slcci  that  was  used 

in  the  wave  propagat  iuu  tests;  iiiterrept  -  a  -  /.«  )  h-.-,/ i  u'' ;  sldpe  -  b  -  i  Jh-s/ln^ 


Damping 


n  =  0.4200 
k  =  0.00373 
r  =  0 
r  =  0 
y^=  0.03611 
g  =  386.4 
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Figure  7.5b. 
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D31F1  —  Incident  Bar 
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7  •  6a .  Input  force  data  t  rum  a  VilPH  test  on  ilry  ptn'ous  sL^iiiilcsb  .steel  Vi/itli  a  10— inch  striker 
bar. 
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Figure  7.6c.  Circumferential  strain  data  from  a  WPB  lest  on  dry  porous  stainless  steel  with  a  10-inch 
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Figure  7.7b. 
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Figure  7.7c. 


D31M1  -  Station  A  (5") 
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Figure  7.9b.  Axial  strain  data  from  a  WPB  test  on  saturated  jjorons  stainJess  steel  wJtJi  a  JO-lnch 
striker  bar. 


00000 


CV2 

CO 

o 

TfO 

o 

O 

o 

o 

O 

o 

o 

o 

o 

o 

o 

Ci 

d 

ci 

I 


^isi:i.U9jajmnojT3 


7-39 


Figure  7.9c. 


D31H1  -  Station 
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Figure  7.10a.  Input  force  data  from  a  WPB  test  on  saturated  porous  stainless  steel  with  a  20-lnch 
striker  bar. 
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Figure  7.10c.  Circumferential  strain  data  from  a  WPB  test  on  saturated  porous  stainless 
a  20-inch  striker  bar. 
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Figure  7.11a.  Input  force  data  from  a  WPB  test  on  satiiruted  porous  stainless  steel  witli  a  20-lnch 
striker  bar  and  a  1/8  Inch  rubber  pad. 
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Figure 


D31R1  -  Station  A  {5") 

D31R1  -  Station  B  (10") 

D31R1  -  Station  C  (15") 
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Figure  7.11c.  Circumferential  strain  diita  from  a  WJ’B  lest  on  saimatuil  porous  stainless  steel 
a  20-incli  striker  bar  ami  a  1/H  imli  rubber  pad. 
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D31M1  -  Saturated 


Figure  7.12a.  Comparison  between 
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-0.000 
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Figure  7.13a.  Comparison  between  axial  strain  records  t  rom  tests  con<hn:ted  i)n  dry  ami  l  luld- 
porous  stainless  steel  specimens  with  a  lO-iiicii  slriUiT  b.ir. 
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Figure  7.13b. 
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Figure  7.14a.  Comparison  between  axial  strain  records  1 roin  tests  conducted  on  dry  and  fluld-sa 
porous  stainless  steel  specimens  with  a  2()-incli  striker  bar. 
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Input  Velocity  Time  History 
MPDAP  Simulation  WFR-V 


Igure  7.15.  Velocity  history  that  was  used  as  the  loading  function  lor  (lie  nninerical  Simula 
test  D31K1. 
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Figure  7.17a. 
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Figure  7.18.  Comparison  of  the  pore  pressure  at  the  5-incli  pt)itU  (Scut  io 
simulation  with  the  corresponding  ipiantity  1 rom  test  UilKI. 


Wave  Propagation  Bar 
MPDAP  Simulation  WFR-\ 
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Propagation  Par.  MPDAP  Siimilation  WPR 
Range  =  5  in,  Mid -rad ins  of  Par 


Figure  7.20a.  Skeleton  velocity  and  apparent  relative  llui.l  velocity  luniputed  by  tlie  miiiierlcal 
simulation  at  tlie  3-iiu:li  point  (Station  A). 
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Figure  7.20c.  Skeleton  velocity  and  apparent  relative  iJuid  velocity  coiiipoted  hy  tlie  numerical 
simulation  at  the  15-lnch  point  (Station  . 
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Position  B. 


Figure  7.21.  Comparison  of  axial  strains  coinpuCcd  by  the  mnnuric 
at  various  radial  distances  from  the  center  of  the 


Wave  Propagation  Bar,  MPDAP  Simulation  WFR-V 
Position  B,  Range  =  10  in 
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